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GENAH USER’S GUTOE AND 
REFERENCE MANUAL 
VERSION 1.0 


Chapter 1 
Introduction 


Genah is the GEneralized Nearfield Acoustical Holography processing sys¬ 
tem. It is used to project acoustical pressure holograms taken in cylindrical 
coordinates to reconstruction surfaces both inside and outside the original 
hologram surface. It provides reconstructions in both the spatial and wave- 
number domains, and its outputs include pressure, normal (radial) velocity 
and intensity. 

First developed in 1985, Genah is a direct outgrowth of planar nearfield 
holographic algorithms whose code dates back to 1980. The code has under¬ 
gone many transitions in its lifetime. It was ported from its original PDP-11 
progenitor under RT-11 to a VAX 11-785 running VMS where it underwent 
extensive modification to (approximately) the form seen today. In 1990, it 
was again ported, this time to a Silicon Graphics platform running Unix. Be¬ 
cause of major differences in both the native number representations and the 
FORTRAN compilers, this latter port engendered a number of major changes 
to the subroutines that support the core Genah algorithm. This allowed 
the code, running on an SGI machine, to read files in DEC number format 
if necessary for backward compatibility. The version included in this 1994 
public release has undergone extensive (although mostly cosmetic) revision 
to provide a uniform user interface and to make the code more readable. 
Nonetheless, the interested user will recognize aspects of coding that harken 
back to the original (very limited) code used on the PDP-11. 

The User’s Guide is intended to provide the user/operator with the infor¬ 
mation necessary to install and intelligently operate the software. It opens 
with a short guide to installing the software on any Unix system with a 

Manuscript approved December 9, 1994. 
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FORTRAN-77 compiler. A very brief tour of Genah is then provided in the 
form of a sample run. Next, a thorough discussion of the primary process¬ 
ing options is undertaken. The user must understand which “direction” the 
processing is to proceed before intelligent parameter choices can be made. 
A more detailed description of the user inputs follows, with information on 
choosing the many parameters necessary to guide Genah ’s processing. Fi¬ 
nally, the input and output files used by Genah are described. This descrip¬ 
tion is meant to allow the user to proceed with operating the software; very 
detailed descriptions of the files, their format and header specifications are 
provided as appendices. 

The Reference Manual is intended to provide the user with an in-depth 
understanding of Genah . The theory of Nearfield Acoustical Holography 
is introduced with references to the literature included. The specialization 
to cylindrical coordinates is demonstrated, including the Green’s Function 
propagator and differences in interpretation of the wavenumber domain from 
planar systems. The implementation of the propagation algorithm is then 
discussed. Finally, a very thorough description is presented of the actual 
code included in this public release. This latter should be sufficient to allow 
for local modification of the code as needed to alter I/O, processing options, 
and so on. 
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Chapter 2 

Installing GEN AH 

2.1 Routines in this Release 

The files in this release consist of a set of FORTRAN-77 routines, a Makefile 
for compiling and linking them, a sample variable filter file, a sample input 
hologram and a sample script file. Table 2.1 lists the files as they appear in 
the distribution (that is, in alphabetical order). Table 2.2 breaks these same 
files down into categories based on their purpose. Those files ending with a 
“. f ” subscript are FORTRAN routines. 

To install Genah , simply make a suitable subdirectory in your filesystem 
and place the files there. If they were saved as one file from an E-mail file, 
you will need to use an editor to separate them. If they were loaded as a 
“tar” (backup tape) file, simply use the “tax -xvf” command to extract 
them. Once a directory listing shows the same files as those listed in Table 
2.1, they are ready to be compiled and tested. 

Genah depends on a set of Bessel function routines to compute the 
propagators. The routines currently used at NRL are those provided by the 
VAXMATH collection of FORTRAN scientific subroutines, themselves descen¬ 
dants of the venerable SLATEC routines. All of the Bessel function routines 
and their support routines are included in this release. They should be com¬ 
piled into an archive (library), referred to in the Makefile as “bessel.a”. (If 
your library will be named diff'erently, e.g. “BESSEL. OLB” , then the make 
file should be altered accordingly.) If your installation includes the double¬ 
precision VAXMATH routines (compiled G-floating on VAXes), then they 
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Bytes 

File Name 

1231 

Makefile 

600 

Sample.Script 

1.4M 

Test.hoi 

46 

Test.window 

905 

CeOlg.f 

1077 

cvt_2byte_to_int4.f 

1187 

cvt_4byte_to_real4.f 

1081 

cvt_decf_ieeef .f 

1131 

cvt_int4_to_2byte. f 

1345 

cvt_real4_to_4byte. f 

2474 

fft.f 

5619 

filread.ieee.f 

3571 

filwrite.ieee.f 

40443 

genah. f 

983 

banning,f 

2446 

header-ieee.f 

4697 

read_f rom_ihead. f 

1443 

tukey.f 

1159 

windo.f 

4677 

write_to_ihead.f 

2467 

wtrev.f 


Table 2.1: List of files in this release of Genah , and their sizes in bytes. 


may be substituted for the Bessel routines provided here. The Makefile (or 
compile command file) should be edited to reflect this before attempting to 
compile and link the files. 

Important Note: The files IlMACH.f and DlMACH.f provided in this 
release are currently set to execute on lEEE.MOST-SIG-BYTE-FIRST-style 
CPUs, such as the SGI MIPS or Sun SPARC. If you will be using some other 
platform, you must edit these two files to reflect the CPU type on which the 
routines will be compiled and executed. 
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Purpose 

File Name 

Main Processing Routine 

genah.f 

Header Include File 

Data File I/O 

header_ieee.f 
filread-ieee.f 
filwrite.ieee.f 

Data File Header 

read_f rom_ihead. f 
write_to_ihead.f 

FP Number Conversion 

cvt J2byte_to_int4. f 
cvt_4byte_to_real4. f 
cvt_int4_to_2byte. f 
cvt_real4_to_4byte. f 
cvt-decf.ieeef.f 

Spatial Windowing 

hainning. f 
tukey.f 

Wavenumber Filtering 

windo.f 

Fourier Transforms 

wtrev.f 

fft.f 

Phase of Complex Number 

cang.f 


Table 2.2: Genah files, grouped by purpose. 

2.2 Compiling Genah : the Makefile 

Genah is composed of FORTRAN routines. It will only execute if compiled 
with a FORTRAN-77 or compatible compiler. It does not use any system- 
specific calls and is self-contained for execution on most systems. 

The compilation of Genah has been greatly simplified through the in¬ 
troduction of a “Makefile.” If the system on which Genah will run has 
the “Make” utility, simply enter make genah at the command line in your 
designated Genah subdirectory. The routines will then be compiled with 
the native FORTRAN-77 compiler and linked (or “loaded”) into an executable 
image called “genah”. On a VAX system, the “FC” variable in the Makefile 
will have to be changed to “FORT”, and the “OPT” variable will have to be 
cleared (no options are needed by the VAX compiler). 

If your system does not have the “Make” utility, you will have to compile 
and link the routines by hand (or with a script file). If you are using a 
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VAX system, no flags are needed. Simply use the FORT command on each 
file, then a LINK command with all the object files listed (with “GENAH.OBJ” 
first) separated by commas. This process can be simplified with a “.COM” 
file, which the user can base on the included Makefile and build using an 
editor. 

If you are using a Unix system, use the f 77 (or its equivalent) command 
with a flag to generate object files (often “-c”) and a flag to force static mem¬ 
ory allocation (often “-static”). This latter flag is not strictly necessary, 
but it has been found to be the safest way to compile the Genah code. The 
object files are then linked (or loaded) with the f 77 command again, using 
the “-o” option to generate the output name “genedi”. Alternatively, some 
systems use the Id command to link object files into an executable image. 
Use whichever is appropriate for your system. Again, see the Makefile for 
syntax and filenames; a shell script could be developed from it. 

Note that Genah has never been compiled above so-called “Level 1” 
optimization (where loop variables are scrutinized and some are placed in 
registers). We cannot guarantee the safe and accurate operation of Genah 
for higher levels of optimization. There is no inherent reason why Genah 
should not perform well under such optimization, but doing so would be the 
responsibility of the user. 


2.3 Testing the Genah Installation 

Once the code has been extracted and compiled, it can be tested with the 
example files included in the distribution. In particular, an example Unix 
script file (Sample.script) is included to run the Genah code, as well as a 
small hologram file (Test. hoi) on which it operates. The user can simply 
execute the script from the command line or use it as a guide to run Genah 
from within a screen shell to get a feel for using it. 

Unfortunately, since hologram files are encoded in a particular binary 
format (see Appendix A), it is not a trivial matter to examine, list, plot or 
visualize them. The user can examine the ASCII output generated for each 
bin in processing Test .hoi as a first-order check on the Genah installation. 
Several plotting and visualization routines have been developed at NRL; they 
are not included in this software release because they are specific to the SGI 
4D family of graphics workstations. Two options are available. First, the 
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user could check with the authors to see if any of the currently available 
tools might be applicable to the user’s system. Second, the user could (and 
probably should) locally develop plotting and listing software, using the I/O 
routines provided in this release (see Appendix D) to load holograms and 
plot or list their results. 
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Chapter 3 

A Quick Tour of GENAH 

The following is a very brief introduction to the Genah software. In Section 
3.1, the structure of the Genah code is outlined in a manner that will famil¬ 
iarize the user with some of the choices to be made for the input parameters. 
In Section 3.2, a sample session with Genah is listed and described. This 
session uses the same inputs as that contained in the script Sample. script, 
provided in this release. 


3.1 Basic Genah Structure 

For detailed information on how the FORTRAN code in Genah is laid out, refer 
to Chapter 9. For information on how that code implements the Genah 
algorithm, refer to Chapter 8. For an in-depth discussion of holography 
theory and how it leads to an algorithm, see Chapter 7. All of these Chapters 
are in the Reference Manual. 

Genah is composed of two primary sections, a preamble and the main 
frequency loop. Within the preamble, the user opens the input hologram 
pressure file and defines all of the flags and parameters needed to execute 
Genah on the input data. The main loop steps through the holograms and 
processes each one according to the directives defined in the preamble. The 
resulting reconstructions are then written to the appropriate output files. 

Opening the input data file sets most of the parameters needed for holo¬ 
gram processing. These include the size of the hologram mesh itself, the 
number of holograms that are available for processing, and the center fre- 
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quencies at which the holograms fall. This information is contained in a file 
“header”, a block of bytes that precedes the data values in the data file. This 
header is fully described in Appendix B. The format of the data files used by 
Genah is described in Appendix A. These parameters serve as guidelines; 
they can be altered before processing proceeds. 

Many decisions need to be made in processing holograms. These options 
are provided as fiags or additional parameters which are set by the user during 
the preamble. They include signal processing capabilities such as window¬ 
ing, filtering, and zero-padding, as well as operational parameters such as 
which holograms to process. The most important parameter is the Process¬ 
ing Option. This option determines the processing path taken by Genah in 
propagating and reconstructing the data. It includes “full” reconstructions 
of pressure, velocity and intensity as well as partial reconstructions in the 
wavenumber domain and in mixed wavenumber-spatial domains. For a com¬ 
plete discussion on the Processing Option, see Chapter 4 in this Guide. For 
more information on the processing directives and their effects on process¬ 
ing, see Chapter 5. All of the inputs and their allowable values are listed in 
Appendix C. 

Once the input file has been opened, its header loaded, and the processing 
directives set, the output files are opened and their headers are written based 
on the input header and the state of the processing directives. The processing 
is then ready to begin. 

The outermost loop of Genah is over individual holograms. Each one 
exists at an independent frequency bin. (Note here the implication for par¬ 
allelization on shared-memory MIMD systems or Parallel Virtual Machines: 
each bin can be processed independently of all other bins.) Within this loop, 
the processing of holograms is carried out in a series of well defined steps. 
This overview will cover all of the steps as if a full reconstruction were re¬ 
quested. Other Processing Options will not use all these steps. 

On each execution of the loop (a new hologram at the next frequency), 
the appropriate hologram is read off the disk from the input file. The first 
operation is the two-dimensional forward Fourier transform of the hologram 
data, moving it from the spatial domain to the wavenumber domain. Since 
the spatial coordinates over which this is done are independent, the transform 
is carried out as the linear combination of one-dimensional transforms, first 
along the axial coordinate and then along the circumferential coordinate. 
The transforms are executed using Radix-2 Fast Fourier Transforms (FFT’s). 
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This has strong implications for the size of the hologram mesh: each direction 
must have a number of points equal to a power of two. Prior to transforming, 
the data can be windowed to ameliorate aperture effects (see Section 8.3). 
The results from the transforms are rotated in the wavenumber domain so 
that the zero wavenumber component resides in the center of the data array. 
This makes for ease of viewing if the processing stops at this step. 

Next, the propagation functions are computed. This involves computing 
ratios of Hankel functions, requiring double-precision computation of vectors 
of Bessel functions. This is one of the more expensive operations in the 
Gen AH algorithm, but there is no good way around it. The functions are 
frequency-dependent, so it doesn’t usually help to compute them offline and 
store them (the disk access would probably exceed the compute time). Once 
computed, they are multiplied with the wavenumber domain representation 
of the hologram, forming an inner product that replaces the expensive spatial 
convolution from which the propagation process is defined (see the theoretical 
discussion in Chapter 7 of the Reference Manual). The data have now been 
“propagated” to a new surface, the core purpose of holographic processing. 
Within the product loop, the data can also be wavenumber-filtered. The 
filter cutoff can be a constant or a function of frequency. In the latter case, 
filter parameters are loaded from an external file before processing begins 
(see Appendix E). There are two propagators, one for pressure and one for 
normal velocity. At the end of this operation, there are now two sets of data 
in the wavenumber domain. 

The data are then ready to be inverse-transformed back to the spatial 
domain. For full reconstructions (both pressure and velocity), two sets of 
transforms are required. The same routines and techniques are used in back- 
transforming as were used in forward-transforming. 

With the data back in the spatial domain, the pressure and velocity can 
be combined via an inner product to form acoustic intensity (see Section 
9.3.11 in the Reference Manual). The total radiated power is the integral of 
the intensity over the reconstruction surface. It is approximated in Genah 
with a rectangular quadrature over the hologram aperture. The intensity 
field is searched for the its maximum value, and all points in the aperture 
that are within 70% of that maximum are listed to an ASCII report file with 
the total power. 

Finally, the reconstructions (pressure, velocity and intensity) can be writ¬ 
ten to their respective output files (see Section 9.3.12). The loop is then 
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repeated for the next hologram of interest. 

The other Processing Options simply step out of this path at one point 
or another and write what could be considered to be intermediate results. 
These can be very valuable, particularly results in the wavenumber domain, 
for other types of analysis. These possibilities are discussed throughout the 
Reference Manual. 


3.2 A Sample Session with Genah 

This section shows an actual interactive session with Genah . It uses the 
inputs contained in the file Sample. script. The user can compare the output 
listings which follow the section with an actual run of Genah for the purposes 
of testing the Genah installation (see Section 2.3). 

First, the full session will be listed out as it would appear on the user’s 
terminal. Then sections of the session will be repeated with brief commen¬ 
tary. For a full description of options and parameters, see Chapter 5 and 
Appendix C. For a complete discussion of processing options and the algo¬ 
rithmic response to parameters, see Chapter 4 and the Reference Manual. 

This session with Genah operates on 45 individual holograms from the 
input data file Test.hoi. A complete normal reconstruction is performed 
using a single wavenumber filter. No additional options are used. 

myhost: Genah> genah 


*** Welcome to GENAH *** 
the Cylindrical Holography Processing Program 
Version 1.0 September 1994 

Input Data Control Switches: 

Conjugate data? (I=exp(iwt)/0=exp(-iwt)): 0 
Enter binary number format (1=DEC/0=IEEE): 0 


Enter the Input Filename: 
Test.hoi 
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Enter the type of units for the input file: 

0 = inches (English) 

1 = meters (MKS) 

2 = centimeters (cgs) 

—> 0 

Enter a choice for axial data window: 

0 = no window (rectaingulax) 

1 = Tukey with 8-point taper 

2 = Hanning (cosine squared) 

3 = Cosine 
—> 1 

Do you wish to ZERO-PAD data axially? 

(Doubles k-space resolution, reduces replicated source error) 
(l=Yes/0=No): 0 

Frequency of first bin = 198.36 

Frequency width of bins = 38.15 

Number of bins available = 45 

Enter the starting hologram #: 1 
Enter the step interval (>=1): 1 
Enter the # of holograms to process: 45 

Enter the k-space filter parameters. 

("rkc" = cutoff wavenumber, "alpha" = Filter slope) 

(Enter " 0, 0" for no filter; 

Enter "-1, 0" to read in a frequency-dependent filter 
from an external file; 

RKC, ALPHA: -1,0 

Enter the file name for the variable filter: 

Test.window 

Additional Circumferential Order Filtering: 
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Enter "1" to remove higher orders, "0" to include all orders: 0 


Select an output processing type: 

0 = Normal Reconstruction (P, V, I) 

1 = K-space velocity withOUT propagation (K) 

2 = K-space velocity WITH propagation (G) 

3 = V(z, n) (NO inverse transform on Phi) (N) 

4 = P(kz, phi) (NO propagation, NO transform on Phi) (Q) 

5 = V(z, phi) for orders "n" = nmin,..., nmax 
—> 0 


For the following target and test: 

3300 Inner Shell Only, Flat End Caps 
TEST HOLOGRAMS FOR GENAH VI.0 (DRHX=0.59055") 
these are the current values in the Header for: 


Axial lattice spacing (inches) = 1.0000 
Shell radius (inches) = 3.2520 
Standoff disteince (inches) = 0.0600 
Hologram Surface Radius (inches) = 3.3120 


Do these values need to be altered? (l=Yes/0=No): 1 
Enter the ajcial lattice spacing, "alat": 1. 

Enter the standoff distance, "drhx": 0.59055 
Enter the reconstruction radius, "radius": 3.25197 


Enter the Output File Name, using the following format: 
"holneime.*xx" where "xx" = 10 x K_r(cutoff), or 
"holname.*vw" where "vw" represents a veiriable filter, 
and the "*" is literal in both cases: 

Test.*vw 
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Bin 

# 

1 

1 

Freq = 

198.4 

1 
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A 

k 

P 

= 
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1 
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163.0 

1 

at 

35 

33 

Bin 

« 

5 

1 

Freq = 

351.0 

1 

Max 

A 

k 

P 

= 

6.1108E+00 

163.7 

1 

at 

38 

34 

Bin 

# 

6 

1 

Freq = 

389.1 

1 

Max 

A 

k 

P 

= 

5.9833E+00 

167.5 

1 

at 

40 

32 

Bin 

# 

7 

1 

Freq = 

427.2 

1 

Max 

A 

k 

P 

= 

6.3682E+00 

-12.1 

1 

at 

38 

51 

Bin 

# 

8 

1 

Freq = 

465.4 

1 

Max 

A 

k 

P 

= 

9.6915E+00 

170.4 

1 

at 

38 

1 

Bin 

# 

9 

1 

Freq = 

503.5 

1 

Max 

A 

k 

P 

= 

1.8090E+01 

176.0 

1 

at 

39 

64 

Bin 

# 

10 

1 

Freq = 

541.7 

1 

Max 

A 

k 

P 

= 

6.2974E+01 

107.5 

1 

at 

25 

1 

Bin 

# 

11 

1 

Freq = 

579.8 

1 

Max 

A 

k 

P 

= 

4.0726E+01 

-20.8 

1 

at 

39 

1 

Bin 

« 

12 

1 

Freq = 

618.0 

1 

Max 

A 

k 

P 

= 

3.2896E+01 

-14.9 

1 

at 

40 

64 

Bin 

# 

13 

1 

Freq = 

656.1 

1 

Max 

A 

k 

P 


5.4091E+01 

161.4 

1 

at 

24 

33 

Bin 

# 

14 

1 

Freq = 

694.3 

1 

Max 

A 

k 

P 

= 

4.8595E+01 

-140.6 

1 

at 

33 

4 

Bin 

# 

15 

1 

Freq = 

732.4 

1 

Max 

A 

k 

P 

= 

2.6263E+01 

162.0 

1 

at 

26 

33 

Bin 

# 

16 

1 

Freq = 

770.6 

1 

Max 

A 

k 

P 

= 

2.1198E+01 

166.2 

1 

at 

24 

33 

Bin 

# 

17 

1 

Freq = 

808.7 

1 

Max 

A 

k 

P 

= 

2.2897E+01 

167.0 

1 

at 

23 

33 

Bin 

« 

18 

1 

Freq = 

846.9 

1 

Max 

A 

k 

P 

= 

5.1654E+01 

-164.1 

1 

at 

22 

33 

Bin 

« 

19 

1 

Freq = 

885.0 

1 

MelX 

A 

k 

P 

= 

2.3355E+01 

-16.2 

1 

at 

33 

1 

Bin 

# 

20 

1 

Freq = 

923.2 

1 

Max 

A 

k 

P 

s 

2.0145E+01 

-10.0 

1 

at 

33 

1 

Bin 

# 

21 

1 

Freq = 

961.3 

1 

Max 

A 

k 

P 

= 

2.1649E+01 

-0.5 

1 

at 

33 

1 

Bin 

« 

22 

1 

Freq = 

999.5 

1 

Max 

A 

k 

P 

= 

2.6507E+01 

19.6 

1 

at 

33 

1 

Bin 

« 

23 

1 

Freq = 

1037.6 

1 

Max 

A 

k 

P 

= 

2.7737E+01 

-39.5 

1 

at 

31 

49 

Bin 

# 

24 

1 

Freq = 

1075.7 

1 

MelX 

A 

k 

P 

= 

3.6167E+01 

-19.9 

1 

at 

30 

17 

Bin 

« 

25 

1 

Freq = 

1113.9 

1 

Max 

A 

k 

P 

= 

1.6432E+02 

167.4 

1 

at 

31 

57 

Bin 

« 

26 

1 

Freq = 

1152.0 

1 

McLX 

A 

k 

P 

= 

1.8496E+02 

178.6 

1 

at 

38 

33 

Bin 

# 

27 

1 

Freq = 

1190.2 

1 

Max 

A 

k 

P 

= 

9.4159E+01 

126.6 

1 

at 

28 

33 

Bin 

# 

28 

1 

Freq = 

1228.3 

1 

Max 

A 

k 

P 

= 

6.4703E+01 

163.6 

1 

at 

21 

1 

Bin 

# 

29 

1 

Freq = 

1266.5 

1 

Max 

A 

k 

P 

= 

5.6077E+01 

-12.5 

1 

at 

44 

1 

Bin 

# 

30 

1 

Freq = 

1304.6 

1 

Max 

A 

k 

P 


4.7581E+01 

-6.5 

1 

at 

44 

1 

Bin 

# 

31 

1 

Freq = 

1342.8 

1 

Max 

A 

k 

P 

= 

5.2260E+01 

4.2 

1 

at 

44 

1 

Bin 

« 

32 

1 

Freq = 

1380.9 

1 

Mcix 

A 

k 

P 

= 

7.7005E+01 

6.7 

1 

at 

44 

1 

Bin 

« 

33 

1 

Freq = 

1419.1 

1 

Max 

A 

k 

P 

= 

8.7289E+01 

163.1 

1 

at 

45 

49 

Bin 

# 

34 

1 

Freq = 

1457.2 

1 

Max 

A 

k 

P 

= 

4.5445E+01 

-177.1 

1 

at 

45 

50 

Bin 

# 

35 

1 

Freq = 

1495.4 

1 

Max 

A 

k 

P 

= 

4.4650E+01 

-162.4 

1 

at 

45 

50 

Bin 

« 

36 

1 

Freq = 

1533.5 

1 

Max 

A 

k 

P 

= 

5.0908E+01 

-141.3 

1 

at 

45 

51 
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Bin 

# 

37 

1 

Freq 

= 

1571.7 

Bin 

# 

38 

1 

Freq 


1609.8 

Bin 

# 

39 

1 

Freq 

= 

1647.9 

Bin 

# 

40 

1 

Freq 

= 

1686.1 

Bin 

# 

41 

1 

Freq 


1724.2 

Bin 

# 

42 

1 

Freq 


1762.4 

Bin 

# 

43 

1 

Freq 

= 

1800.5 

Bin 

« 

44 

1 

Freq 


1838.7 

Bin 

# 

45 

1 

Freq 

= 

1876.8 


myhost: Genah> 


Max A & P = 
Max A & P = 
Max A & P = 
Max A & P = 
Max A & P = 
Max A & P = 
Max A & P = 
Max A & P = 
Max A & P = 


We can now take a closer look at this session. 


6.0880E+01 

92.5 

1 

at 

27 

34 

6.4391E+01 

148.1 

1 

at 

45 

2 

7.4035E+01 

148.2 

1 

at 

26 

33 

5.1020E+01 

156.6 

1 

at 

26 

33 

4.1062E+01 

157.8 

1 

at 

26 

33 

3.8541E+01 

162.0 

1 

at 

25 

33 

3.9626E+01 

165.9 

1 

at 

29 

27 

5.9036E+01 

169.3 

1 

at 

29 

27 

2.0275E+02 

175.3 

1 

at 

29 

27 


Before the file can be 


opened, the software needs to know what sense of time was used in creating 
the hologram reference file and in which number format the data is stored. 
Then the filename can be entered and the file opened. The power calculations 
at the end of Genah require knowledge of the engineering units in which the 
data (pressure per unit force) are stored, so this too is entered: 


Input Data Control Switches: 

Conjugate data? (I=exp(iwt)/0=exp(-iwt)): 0 
Enter binary number format (1=DEC/0=IEEE): 0 


Enter the Input Filename: 
./Test.hoi 


Enter the type of units for the input file: 

0 = inches (English) 

1 = meters (MKS) 

2 = centimeters (cgs) 

—> 0 

A series of questions establishes guidelines for processing the data. These 
include zero-padding the input data before processing, the actual holograms 
to be processed, the wavenumber filter to be used, and whether or not re¬ 
constructions should be restricted to a subset of available circumferential 
harmonics. Note that a frequency-dependent filter is invoked. 


Do you wish to ZERO-PAD data axially? 
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(Doubles k-space resolution, reduces replicated source error) 
(l=Yes/0=No): 0 

Frequency of first bin = 198.36 

Frequency width of bins = 38.15 

Number of bins available = 45 

Enter the starting hologram #; 1 
Enter the step interval (>=1): 1 
Enter the # of holograms to process: 45 

Enter the k-space filter parameters. 

("rkc" = cutoff wavenumber, "alpha" = Filter slope) 

(Enter " 0, 0" for no filter; 

Enter "-1, 0" to read in a frequency-dependent filter 
from aji external file; 

RKC, ALPHA; -1,0 

Enter the file name for the variable filter: 

Test.window 

Additional Circumferential Order Filtering: 

Enter "1" to remove higher orders, "0" to include all orders: 0 

The user must choose the processing path (and consequently, the type 
of output) used by Genah . Here, a full reconstruction, including pressure, 
normal velocity and intensity, is requested. Then the user is asked to review 
the critical input parameters alat (the axial lattice spacing), radius (the 
shell radius), and drhx (the radial distance to the reconstruction surface 
from the hologram surface, or “standoff”). Since the standoff parameter 
drhx appears to have been corrupted in the header, it is corrected here. 
Then a reconstruction from the hologram surface to the shell’s surface can 
be completed. 

Select an output processing type; 

0 = Normal Reconstruction (P, V, I) 

1 = K-space velocity withOUT propagation (K) 

2 = K-space velocity WITH propagation (G) 

3 = V(z, n) (NO inverse transform on Phi) (N) 
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4 = P(kz, phi) (NO propagation, NO transform on Phi) (Q) 

5 = V(z, phi) for orders "n" = nmin,..., nmax 
—> 0 

For the following target and test: 

3300 Inner Shell Only, Flat End Caps 
TEST HOLOGRAMS FOR GENAH VI.0 (DRHX=0.59055") 
these are the current values in the Header for: 


Axial lattice spacing (inches) = 1.0000 
Shell radius (inches) = 3.2520 
Standoff distance (inches) = 0.0600 
Hologram Surface Radius (inches) = 3.3120 


Do these values need to be altered? (l=Yes/0=No): 1 
Enter the axial lattice spacing, "alat": 1. 

Enter the standoff distance, "drhx": 0.59055 
Enter the reconstruction radius, "radius": 3.25197 


Finally, the user must enter an output filename template^ which the soft¬ 
ware uses to open the necessary output files. Notice the use of the extension 
. *vw. This indicates the use of a frequency-dependent filter, while the * is 
used as a marker to substitute the correct extension modifier for each output 
file type. 

Enter the Output File Name, using the following format: 
"holname.’t'xx" where "xx" = 10 x K_r(cutoff), or 
"holname.*vw" where "vw" represents a variable filter, 
and the is literal in both cases: 

Test.*vw 

At this point, processing begins. As each hologram is loaded, it is scanned 
for its maximum modulus, the phase at that point, and the coordinates in 
the lattice where this occurs. These are typed to the screen; they are also 
written to the ASCII record file (with the “.Ovw” filename extension). They 
are not repeated here; refer to the output listing at the beginning of this 
section. 

When finished, Genah will have written four files for this processing 
option. Test. Ovw will contain an ASCII description of the run, with the input 
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file maxima and output file statistics related to the maximum intensity and 
power in each hologram. Test.Pvw and Test.Vvw will contain the pressure 
and normal (radial) velocity reconstructions, respectively, in the hologram 
binary format. Test. Ivw will contain the intensity reconstructions in the 
hologram binary format. These files are now ready to be used in visualization 
and further analysis. 
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Chapter 4 

Processing Options 


4.1 The Gen AH Process Path 

One of the powerful aspects of Genah is that it is built up from a set of 
linear operators. Processing holographic data with these operators can be 
likened to traversing a path along the building blocks of the algorithm, as 
shown in Figure 4.1. Likewise, since each dimension of holographic data is 
independent (two spatial and one time), the Fourier transform operations on 
them can be mixed in different orders. More importantly, the result of any 
particular traversal of the operations path has an interpretation which can 
bring physical information about the system to light. 

The choice of which path to traverse in Genah is, therefore, the most 
critical one in processing holographic data. The sections of this chapter 
are concerned with an introduction to the paths currently available in the 
Genah code. For more detailed information about the nature and impact 
of these paths, see the Reference Manual. It is important that the user 
understand that other path choices could be coded. The ones currently 
available have proven to be the most used at NRL to date; the first three 
choices constitute 99% of all the Genah processing done at NRL. Also, rather 
than pile every conceivable option into the Genah code, other independent 
programs have been developed to take results from Genah and perform other 
operations on them. In one important example, data are returned to the time 
domain from the frequency data in Genah output to make animations of 
impulse responses. Similarly, the conversion of digitized time histories into 
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Operation 

Option 0 

Option 1 

Option 2 

Option 3 

Option 4 

Option 5 

FwdFT(z) 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Fwd FT((j)) 

Yes 

Yes 

Yes 

Yes 

No 

Yes 

Propagate P 

Yes 

No 

Yes 

Yes 

No 

Yes 

Propagate V 

Yes 

No 

Yes 

Yes 

No 

Yes 

Inv FT(k,) P 

Yes 

No 

No 

Yes 

No 

No 

Inv FT(k,) V 

Yes 

No 

No 

Yes 

No 

Yes 

Inv FT(n) P 

Yes 

No 

No 

No 

No 

No 

Inv FT(n) V 

Yes 

No 

No 

No 

No 

Yes* 

Result 

pXHzA) 

p{k,,n) 

Vik,,n) 

V(z,n) 


V{z,4>') 


Table 4.1: Genah Processing Options Summary (* limited n) 


frequency holograms is carried out prior to running data through Genah , 
though it could conceivably be done within Genah . In any case, the options 
currently available in Genah are by no means exhaustive and could certainly 
be augmented. 

The processing options are summarized in Table 4.1. 


4.2 Option 0: Normal Reconstruction 

This is the “workhorse” option in Genah ; it represents 90% of the use of 
Genah code. In “Normal Reconstruction”, the full processing path of Fig¬ 
ure 4.1 is traversed. Pressure holograms are transformed to the wavenumber 
domain, filtered, propagated to a new surface (and converted to velocity) and 
inverse transformed to spatial variables. The result is a set of pressure, nor¬ 
mal velocity, and acoustic intensity holographic reconstructions at a surface 
concentric to, but different from, the original. With the intensity available, 
the total power in the aperture is computed. Locations in the spatial mesh 
which are within a specified percentage (currently 70% ) of the maximum in¬ 
tensity are highlighted, and the relative phases of the pressures and velocities 
at these points give indications of how reactive the sound field is locally. 

This process can involve the use of a spatial tapering window in the axial 
direction to mitigate truncation effects in the data aperture. It can require 
zero-padding in the axial direction, either to pad out to the nearest power of 
two number of points (for the Radix-2 Fast Fourier Transform routine used 
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in Genah ) or to reduce aliasing from the virtual replicated sources that 
arise from the use of a Discrete Fourier Transform. It almost always requires 
the application of a wavenumber filter to reduce the high wavenumber noise 
inherent in backpropagation: this filter may need to change as a function 
of frequency. All of these are presented as choices by the software when a 
normal reconstruction is requested. 

The key parameter in reconstruction is the radial distance between the 
original hologram surface and the reconstruction surface. Since Genah has 
always been used to project data onto cylindrical shells, this distance is usu¬ 
ally set as the “standoff” between the data acquisition probe (hydrophone) 
and the outer radius of the shell. This parameter must be set properly and 
checked in order to correctly reconstruct holograms backward to the shell. 


4.3 Option 1: p(kz^n) without propagation 

As mentioned above in Section 4.2, backpropagated reconstructions almost 
always require wavenumber filtering. The question arises, “How do I choose 
the filter parameters?” Most importantly, for the circular (2-dimensional) 
filter provided in Genah , where should the cutoffs be placed? The quickest 
way to answer this question is run Genah on a small, representative set 
of holograms using this Processing Option 1. This option simply forward 
transforms the data to the wavenumber domain (or k-space, as it is com¬ 
monly known) without propagating to another surface. When the results 
are visualized (mosaics of false-color plots are most useful), the wavenumber 
regions containing meaningful information are usually clearly separated from 
higher wavenumber noise components. The gap between the two is a good 
place to locate the filter cutoff. 

Unpropagated k-space reconstructions are also good at diagnosing any 
potential problems in the data. These might include aliasing, which then 
requires zero-padding, and trunctation effects, indicating the need for pre¬ 
windowing the data spatially in the axial direction. When plotted in the 
frequency-wavenumber domain as a function of circumferential order, dis¬ 
persion surfaces reveal which wavetypes exist in the nearfield of the shell. 
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4.4 Option 2: V{kz^ n) with propagation 

When a filter has been selected, it is best tested in the wavenumber domain. 
This is accomplished using this processing option, which, in addition to the 
forward transform of the previous section, applies the velocity propagator. 
The results, when visualized, will clearly reveal whether the wavenumber 
filter cutoff(s) have been properly set. In this way, this processing option 
provides another diagnostic step along the road to full reconstruction. 

However, it is also extremely important in and of itself in generating 
frequency-wavenumber plots (dispersion surfaces), as mentioned above. Since 
these data have been propagated (usually to the shell surface), they represent 
the actual normal velocity wavenumber content of the shell itself, a key com¬ 
ponent in diagnosing shell physics. Thus, output from this processing option 
is critically important to many analyses of shell vibrations. Further wave- 
number filtering can be developed to limit the data to left- or right-traveling 
waves only or to exclude subsonic or supersonic waves. The results of this 
Processing Option 2 are critically important in assessing the application of 
further operations. 


4.5 Option 3: V{z,n) (no inverse transform 
on (p) 

This is the first of the specialty options in the Genah processing path. Af¬ 
ter having been converted to velocity and propagated to the new surface, 
the data are inverse transformed only in the axial direction. This results in 
reconstructions which are still separated as a function of circumferential or¬ 
der. In effect, the routine just bypasses the transform back to angle. This is 
valuable in examining wave activity specific to certain orders. Many scatter¬ 
ing and radiation phenomena are likewise classified by these circumferential 
orders. The lowest orders (0, 1 and 2) are often of greater importance in 
the study of shell vibrations, since they contribute most of the power to the 
acoustic farfield. 

Processing Option 3 is also valuable in comparing with the results of many 
available Finite/Infinite Element codes for axisymmetric bodies. These codes 
divide their processing across the independent orders as well. In some cases. 
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it is more important to compare the synthetic models to the experimental 
data order by order rather than having them summed. Option 3 enables 
these comparisons. 


4.6 Option 4: p{kz^ (/>) (no propagation, no 
transform on 0) 

This is the second of the specialty options in the Genah processing path. 

It is not one of the most useful, but it is definitely the simplest. The result 
of this option is simply to forward transform the pressure holograms in the 
axial direction only. The results are then written as a function of angle at the 
original hologram surface, that is, No propagation or wavenumber 

filtering is performed. 

Option 4 can be used on non-holographic data which require one-dimensional 
forward Fourier Transforms. For instance, if such data were stored in a single 
“hologram” as a function of one spatial dimension and frequency, this option 
could be used to transform the spatial dimension to wavenumber without 
modifying the data in any other way. Thus, Genah can be used as a tool 
for transforming data that is not actually holographic. 


4.7 Option 5: 14 ( 2 :, (^'), Limited Orders 

This is the third of the specialty options in the Genah processing path. It is 
the most similar to a normal reconstruction, with two exceptions. First, only 
velocity is reconstructed; pressure and intensity are not available. Second, 
the inverse transform from circumferential orders to angle is limited to a 
subset of the available orders as specified by the user. Unlike the filtering 
option which limits the orders to (—n„,aj.,...,+n,nax;Oinclusive), Option 5 
allows the user to specify a range of orders n = {nmini • • • > ^max) inclusive. 
Thus, the user could effectively exclude the lowest circumferential order(s) if 
so desired. 

However, the remaining orders axe not Fourier summed into a single re¬ 
construction at each frequency. Instead, for each order |n|, both the positive 
and negative components are summed back to an angular representation in 
(j). The reconstruction at each order |n| is then written for each frequency. 
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Thus, if M bins are processed with N orders each, the resulting file will 
contain M x N reconstructions. This will impact any post-processing rou¬ 
tines which compute a given reconstruction’s frequency by multiplying the 
bin number by the frequency interval. 
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Figure 4.1: Schematic of the Nearfield Acoustical Holography processing 
path. 
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Chapter 5 

Inputs and Parameters 


A very brief introduction to the user inputs needed for a Genah run was 
given in Chapter 3. This chapter is intended to explicate all of the possible 
responses to requests for inputs in Genah and how these responses will 
affect the outcome of the processing. It is not intended as a tutorial on the 
algorithm: that information is included in Chapter 8. Instead, this chapter 
will simply define allowed responses and indicate how the routine will respond 
to the choices the user makes. 


5.1 Input Data Information 

Each file written in NAH format (see Appendix A) has a short header (see 
Appendix B) with a good deal of information about the data in the file. 
Nonetheless, some information must be known about the file before it can be 
located, opened and read. This section describes four items (and one error 
check) that Genah needs to accomplish this. 

5.1.1 Input Data Control Switches 

The following two data type switches are included in the current version of 
Genah for historical reasons. 

The first switch (iconj) allows the user to conjugate all incoming data. 
This has the effect of reversing the sense of time in the expression of a prop¬ 
agating wave, that is, either or The Fourier transforms 
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within Genah are fixed to recognize the first sense as an outgoing (radially) 
or rightward going (axially) wave. However, different preprocessing routines 
which generate the input holograms from digitized data may use different 
FFT routines with confiicting senses of time (or in this case, forward vs. in¬ 
verse transforms). In particular, the FFT routines provided with MATLAB 
are known to cause this problem. The user should find out which sense of 
the imaginary number is used in any preprocessing. If this is not possible, it 
can be tested by processing a few higher frequency holograms and examining 
data well outside the test region for increasing phase in the near to far field 
transition (outward going waves). If the phase of the reconstruction is gener¬ 
ally decreasing away from the source, then the opposite time sense must be 
employed. When the FFT’s provided with Genah are used in preprocessing 
(and the time and spatial transforms use opposite senses, respectively), then 
respond to the first switch with “0”. Otherwise, respond with “1”. Again, 
be sure to test any new situations to be sure. 

Summary: iconj Conjugate data?(l=exp(iwt)/0=exp(-iwt)): 

Response: 1 Conjugate all incoming data. 

0 Do not conjugate any incoming data. 

When running on a machine which uses the IEEE number format (most 
workstations and processors). The second switch (ivax) allows Genah to 
read in binary hologram data stored in DEC number format. It simply con¬ 
verts all incoming values under the assumption that the data are DEC binary 
while the executable is running under IEEE binary. Note that the converse 
is not provided as an option in this release of Genah (ie., conversion from 
IEEE format to DEC format). If such a need were to arise, it would require 
modification of the source code. Contact the authors in that case. This 
switch was provided because the code, originally running on VAX machines, 
was ported to SGI machines, while data taken up to that time was not. Note 
also, that if both the data and the code are native to a DEC machine, the 
user would respond as if no conversion was necessary (i.e., using the IEEE 
switch). Note that some machines, such as the Convex C-Series, have both 
native (IEEE) and DEC-compatibility modes. In such cases, it is advisable 
to compile and run Genah in the mode in which data is stored rather than 
using conversions. 

Summary: ivax Enter binary number format (1=DEC/0=IEEE): 
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Response: 1 
0 


Convert from DEC format to IEEE format. 

Do not perform any number format conversions. 


5.1.2 Input Data File Name 

The name of the disk file containing holographic input data (f name) includes 
the full directory path and the file’s native name. Note that the string as 
entered must not exceed 80 characters in length. This can easily occur when 
data is being read from a remote disk and/or from deep within a directory 
tree. Care should be exercised in these cases. When a conflict arises, it is 
best to run Genah from a directory nearer to the data (to shorten the name) 
or to set up an alias (Unix) or symbol (VMS) which points to the data. 


Summary: fname 
Response: (Unix) 


(VMS) 


Enter the Input Filename: 

name or 

./name or 

dir/name or 

dir/subdir/.. ./name or 

/dir/subdir/.. ./name or 

/partition/dir/subdir/.. ./name etc. 

name.ext or 

[dir]name.ext or 

[dir.subdir...]name.ext or 

disk: [dir.subdir...]name.ext or 

machine:: disk: [dir. subdir...] name. ext etc. 


5.1.3 Engineering units for input spatial measures 

Data coming into Genah is assumed to be acoustic pressure per unit ref¬ 
erence (a force gauge, an accelerometer, another hydrophone or the drive 
voltage) in SI units, e.g. Pascals per Newton. However, the engineering 
units used for spatial measure of the aperture and standoff can be in any 
system of units, because the propagator functions depend only on the ratio 
of spatial measures. (An exception to this is the specific acoustic impedance 
pc, used in converting from pressure to velocity. It is fixed in SI units for 
water at atmospheric pressure at 20°C.) Even so, the units of spatial measure 
must be known to the routine for two reasons. First, wavenumber filtering 
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requires knowledge of the axial lattice spacing to form axial wavenumbers. 
This impacts the values entered in for the wavenumber filter parameter rkc 
{K,cutoff) as well. Second, computation of total power requires integration 
of the intensity over the total aperture, whose area must be known. 

Historically, the lattice spacing and standoff were measured exclusively in 
inches. More recently, these values have been stored in meters or centimeters 
for compatibility with other acoustic measurement programs. Since an entry 
was not provided in the file header (see Appendix B) for a units flag, this 
has been left as an input parameter. 

Summary: units Enter the type of units for the input file: 

Response: 0 inches (English) 

1 meters (MKS) 

2 centimeters (cgs) 

5.1.4 Correction for number of points in the circum¬ 
ferential direction 

When NAH data are taken at NRL, an extra axial scan line is always added 
which overlays the first scan line. In this way, data are available from the 
very beginning and end of the experiment which should theoretically be iden¬ 
tical. This provides a good measure of the variation which has occured over 
the course of the experiment. This variation commonly arises from water 
temperature changes or amplifier level drifting. The number of axial scan 
lines surrounding the source, i.e. the number of lattice points in the circum¬ 
ferential direction (ncirpts), is always chosen to be even and almost always 
a power of two for the radix-2 FFT in Genah . The extra line is counted in 
before preprocessing. It can appear in the original data header (see B as an 
odd number of points in the circumferential direction, and there have been 
instances in which this is not corrected in generating the input hologram file. 

Since Genah depends on this number for all of its operations, especially 
the FFTs, it must be even. (The only exception is in the event of a single scan 
line, in which case the data are replicated to form a synthetic, axisymmetric 
field.) In the event that this value is not an even number when Genah opens 
the header field for the input file, an error is flagged, and the user is given 
the opportunity to correct the value. Note that the actual number of data 
lines in each hologram must agree with the corrected value; otherwise, the 
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algorithm for counting records in the file will fail. Again, this input will only 
be requested if Genah detects an odd number in the header location. The 
user must know if this situation will occur before writing a shell script to 
execute Genah or the subsequent inputs will be out of order and cause the 
routine to fail. 

Summary: ncirpts *** ERROR *** in "ncirpts" 

Current value is NOT acceptable. 

Please reenter "ncirpts": 

Response: “n” (even) Actual number of scan lines in file 

(must be even, power of 2) 


5.2 Processing directives 

This section describes seven parameters which control the processing done 
by Genah . These include spatial windowing and padding, choosing which 
holograms (frequency bins) to process, wavenumber-space filtering, process¬ 
ing path option and reconstruction parameters. 

5.2.1 Axial direction spatial window 

Hologram measurement apertures axe usually designed to extend sufficiently 
past the ends of the source that the acoustic pressure at the edge of the 
aperture has dropped below the measurement signal to noise ratio (SNR). 
This is referred to as “overscan”. For more details regarding this, see Section 

8.3 in the Reference Manual. For many measurements, this cannot quite be 
achieved, either because full overscan would not allow a sufficient number 
of points to fall over the source itself, or because at certain frequencies the 
source tends to project much more energy axially than over the rest of the 
frequency range of interest. For whatever reason, any significant pressure 
level remaining at the ends of the aperture will appear as a truncation dis¬ 
continuity to the Fourier Transforms. This can strongly bias the results of a 
reconstruction. To avoid this, four tapering windows (one of which is simply 
the rectangular window) have been made available in Genah to force the 
data to zero at the ends of the aperture. 

The first choice is simply not to taper the data at all. This rectangular 
or “box-car” window is fine if the overscan is sufficient at all frequencies. 
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It yields the narrowest main lobes possible in the wavenumber domain, but 
with the highest sidelobes. 

The second choice is the workhorse for most NRL processing. This is the 
Tukey window. It has no taper over most of the aperture (that is, a value of 
1.0), then a short cosine taper over the last eight points at each end. This 
yields wavenumber responses almost identical to those of the rectangular 
window for apertures of 128 or more points while eliminating truncation 
problems. 

The third and fourth choices are cosine tapers. The Hanning, or cosine- 
squared taper is the most severe of all the windows. The Hamming, or 
cosine window is less severe. Both windows affect the entire aperture, with 
a value of 1.0 only at the center point. The Hanning window provides one 
more continuous derivative at the boundary than the Hamming, yielding 
faster sidelobe roll-off in the wavenumber domain. These windows must be 
used with caution. They tend to bias the reconstructions strongly (they leave 
their footprint, so to speak, on the results). They are most useful in reducing 
background in the wavenumber domain; this is critical in trying to find weak 
trends in visualizations of that domain. Their use at NRL has been confined 
to sources with very strong endfire and in scattering holography with bow 
incident wave energy. 

Note that adding windows of other types to the Genah code would be 
straightforward. The user would simply add it to the list of choices and 
place the appropriate subroutine call into the if ... then tree which loads 
the windowing vector, tuk. 


Summary: 

Response; 


window-choice Enter a choice for axial data window: 
0 No window (Rectangular) 

1 Tukey with 8-point taper 

2 Hanning (cosine squared) 

3 Hamming (cosine) 


5.2.2 Axial direction zero-padding 

Much of the art of processing nearfield holograms lies in taming the artifacts 
introduced by the use of discrete, finite length Fourier transforms. One 
powerful weapon in the user’s arsenal is the ability to zero-pad the input 
stream. This has two valuable effects. First, it doubles the distance between 
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the true source and its replicated images, greatly reducing the effects of 
overlap aliasing (i.e., the energy from the image sources is less likely to leak 
into the source aperture). Second, padding in the spatial domain has the 
effect of interpolating in the wavenumber domain. While no true additional 
information is introduced (the interpolated points fall on sin{x)/x curves 
whose peaks occur at the original wavenumber values), the interpolation is 
very helpful in making plots of wavenumber response. 

For these reasons, Genah includes an option to zero-pad the input data. 
The option should be used with care, as it increases processing time and 
memory usage. Output holograms in the wavenumber domains (Processing 
Options 1 and 2) will be twice as large as the input holograms. Note that 
reconstructions back to the spatial domain do not retain the extra points; 
they are written out using the same number of axial lattice points as the 
input files. 

Summary: izadd Do you wish to ZERO-PAD data axially? 

Response: 1 Yes, zero-pad input data 

0 No, do not zero-pad input data 


5.2.3 Hologram frequency bin control 

It is not always the case that the user wishes to process every available 
hologram in an input file. A set of controls is made available to select the 
sequence of a subset of the holograms. These include the bin number at 
which the sequence begins, the step interval (i.e., the number of bins to 
skip between holograms, including the next to be processed), and the total 
number of bins to process in this run. It is worthwhile to insure familiarity 
with this method of enumerating which holograms are to be processed from 
the original set. An example case will make this more clear. 

Suppose 1000 holograms are available in the original data set. The first 
hologram is associated with the frequency bin centered at 1000 Hz, and each 
bin has a width of 10 Hz. This means the last available bin will be centered 
at 10990 Hz. If the user wishes to process all of these holograms, the starting 
hologram will be “1” (even though it is not at the O-frequency bin), the 
step interval will be “1” (incrementing to each successive bin) and the total 
number of holograms to process would be “1000”. If, however, the user needs 
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information only on the range from 2000 to 3000 Hz in 100 Hz increments, 
the starting hologram would be “101”, the step interval would be “10” and 
the number of holograms to process would be “101”. 

The most important distinction here is that one does not enter the first 
and last holograms to process, but the total number of holograms to be 
processed. Also of importance is the fact that if any processing were to be 
done on the resulting (processed) data, the bin at 2000 Hz would become the 
first bin of that data, and a step interval of 1 would span 100 Hz. In other 
words, hologram counting is always local to the file being processed; there is 
no “memory” of which bins were selected from a previous input file. 


Summary: 

Response: 

ihStcLTt 

“f’ 

Enter the starting hologram #: 

The number of the first frequency bin to be 
processed from the available set 

Summary: 

Response: 

ihstep 

“m” 

Enter the step interval (>=1): 

The increment to the next successive frequency bin 
to be processed from the available set 

Summary: 

Response: 

ihno 

“n” 

Enter the # of holograms to process: 

The total number of frequency bins to be 
processed from the available set 


5.2.4 Wavenumber filter control parameters 

Wavenumber filtering is a key task in nearfield acoustical holography pro¬ 
cessing. It is perhaps the one aspect of the processing that continues to 
require human intervention, eluding fully automated or “turn-key” opera¬ 
tion. Of all the parameters available to the holographer, it is the filters that 
most infiuence the outcome of reconstruction. The need for filtering and the 
theoretical underpinning for it are more fully explored in Chapter 7 of the 
Reference Manual. The process of actually selecting the cutoff and slope 
parameters are described in Appendix E. This section confines itself to the 
mechanics of entering the filter parameters into the Genah algorithm. 
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Only a single type of filter is available in the current implementation 
of Genah . It is a fiat-top filter with a matched exponential taper at the 
edge. Researchers have explored various other filters, but no particular ben¬ 
efits were established. The current filter is well understood, robust and has 
withstood the test of time. 

The filter is computed live within Genah , and at any given frequency 
only two parameters are needed. The first, kcutoff-, is the wavenumber at 
which the taper, or “cutofi’”, begins. This value is measured in radians 
per unit length, or wavenumber. The axial wavenumber is computed as 

= 27r/A^, and the circumferential wavenumber is k^ = n/r. (Note that 
while is continuous, k^ exists only at the integer values of n. However, k^ 
becomes discretized in practice when the aperture is sampled. See Chapter 
8 for more information.) This requires that the user know what units are 
being used to measure axial lattice point spacing and the radius (see Section 
5.1.3 above). The second parameter, a, controls the slope, or rate of taper, 
of the exponential tail. This parameter will be explored further in Appendix 
E; for now, it is sufficient to state that a « 0.05 is acceptable. 

Since the region of useful wavenumbers tends to grow with frequency, a 
single filter setting may not be adequate for the entire range of frequencies to 
be processed. In this case, an option is available to read in a set of filters from 
an external ASCII file. Again, for details, consult Appendix E. Also, in order 
to explore where the cutoflf should be placed, a preliminary reconstruction of 
velocity wavenumbers with no filtering may be needed. An option is available 
to apply no filter at all. 


Summary: rkc, alpha 
Response: “0, 0” 

“- 1 , 0 ” 

“a, 6” 


Enter the k-space filter parameters 
Do not apply any wavenumber filter 
Read in a frequency-dependent filter from 
an external file 

kr{cutoff) = o radians per length. Slope = b 


If the user chooses to read in a frequency independent filter set from an 
external file, Genah immediately prompts for the name of that file. 


Summary: wname Enter the file name for the variable filter 
Response: “name. ext” The full path and file name of the file which 

contains the parameters for the filter set. 

(Limited to 80 characters or less.) 
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There is not much error checking on the filter parameters. If, however, a 
filter set has more filter points specified than can be stored in the available 
array, the following error will be announced: 

*** FATAL ERROR *** in Variable Filter: 

Number of entries exceeds size of "wbin". 

Increase parameter "nrkcbins" eind recompile to continue. 

The program will be terminated immediately. The user should then edit the 
parameter “nrkcbins” in the Genah code, increasing it enough to accomo¬ 
date the larger filter. In this release of Genah , nrkcbins has been set to a 
value of 250. 

For cylindrical holograms, the circumferential wavenumbers, or orders, 
are independent of one another. They may be additionally filtered (removed 
completely) without impacting other orders or their individual contributions 
to the reconstructions. An extra option is provided in Genah to remove 
(zero) any orders above a chosen maximum, in addition to any filtering in the 
wavenumber domain. Order filtering is available for all Processing Options 
except Option 5 (see Chapter 4), for which separate order filtering is available. 

Note that this general order filtering is low-pass only, i.e., all orders from 
—rimax to +^max (including 0) are automatically included. 

Summary: nfilt Additional Circumferential Order Filtering: 

Response: 0 Include all available circumferential orders 

(No filtering) 

1 Remove circumferential orders above a set 

maximum 

If a response of “1” is given, then Genah further requests the (integer) 
maximum circumferential order to be used in reconstructions. This number 
must be less than or equal to one half the total number of circumferential 
lattice points (or axial scan lines). For example, if there are 64 circumfer¬ 
ential data locations, then the available orders would include —32,..., -1-31, 
inclusive. It would not make sense to request orders higher than 32 to be 
removed in this circumstance. 

Enter the maximum order "n" to INCLUDESummary: nf_max 
Response: “rimax" The maximum circumferential order included 

in reconstruction (< ncirpts/2). 
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5.2.5 Genah processing path option 

The Genah processing path is the key option in choosing a reconstruction 
result type. The nature of these paths was discussed in detail in Chapter 4; 
those details will not be revisited here. This section is included for complete¬ 
ness and as a reference aid for the user. In the table that follows, the capital 
letter following each entry is the identifier added to the filename extension 
when the output file is opened. 


Summary: ifqout 
Response: 0 
1 
2 

3 

4 

5 


Select an output processing type: 

Normal Reconstruction (P, V, I) 

K-space velocity without propagation (K) 

K-space velocity with propagation (G) 

V{z,n) (NO inverse transform on 0) (N) 

P{k 2 , 4 >) (no propagation, no transform on 4>) (Q) 
V (z, 4>') for orders n = ..., Umax- 


If Processing Option 5 is chosen, this is the point at which Genah will 
request the actual range of circumferential orders to be included. Unlike 
the general circumferential order filtering mentioned in the previous Section 
which is low-pass only, both a minimum and a maximum may be specified 
with Option 5. However, only the reconstructed velocity will be available as 
an output. Note that the specified minimum and maximum apply to both 
positive and negative orders. For example, in a scan with 64 circumferential 
lattice points, orders would include —32,..., 4-31, inclusive. If the Option 5 
filter was set at “1,3”, then the actual orders used in reconstruction would 
include (±1, ±2, ±3). This would remove the n = 0 “breathing mode” while 
limiting reconstruction to the (usually more important) modes below the 
fourth. 

Summary: main, nmax Enter the minimum and maximum circumferential 

orders to keep in processing: 

Response: Lowest and highest orders to retain 

(0 < n,nin,n^ax < ncirpts/2.) 


5.2.6 Processing pzirameter adjustment 

There are several key spatial measurements without which holographic recon¬ 
structions cannot proceed. The parameters include the axial lattice spacing. 
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the shell radius and the standoff distance from the shell radius to the holo¬ 
gram measurement surface. These parameters are carried in the hologram 
file header (see Appendix B). However, from time to time they can become 
corrupted prior to processing. They are of sufficient importance that Genah 
always lists them for the user prior to proceeding with reconstructions, and 
ofi'ers the user the opportunity to alter or correct them. 

One particular case in which they might be altered even though they are 
correct is when the user wishes to reconstruct to some surface other than 
the outer radius of the shell. A typical example would be when the normal 
velocity and intensity are needed at the hologram surface itself (a propagation 
distance of zero). In that case, the shell radius would be increased to be equal 
to the hologram surface radius and the standofi" distance would be reduced 
to zero. 

Note that the hologram surface radius is a derived quantity, computed by 
adding the shell radius and the standofi" distance. This is due to the fact that 
the latter two are measurable quantities. Nonetheless, for reconstruction pur¬ 
poses, the hologram surface radius is the radius from which reconstructions 
proceed. If the reconstruction radius needs to be altered from the shell sur¬ 
face, be sure that the sum of that new radius and the new standoff distance 
are equal to the original hologram surface radius. Otherwise, the results will 
be biased. As an example, suppose a 5 inch radius shell is measured with a 
1 inch standoff. The parameters would be announced as follows: 

these axe the current values in the Header for: 


Axial lattice spacing (inches) = 1.0000 
Shell radius (inches) = 5.0000 
Standoff distance (inches) = 1.0000 
Hologram Surface Radius (inches) = 6.0000 


For these conditions, the hologram surface radius is reported as 6 inches, the 
sum of the shell radius and the standoff distance. Suppose the user wished 
to reconstruct only halfway in, say to the surface of a half-inch thick test 
tile. Then the shell radius would be increased to 5.5 inches and the standoff 
distance would be reduced to 0.5 inches. The hologram surface radius would 
remain at 6.0 inches, but Genah would only reconstruct inward 0.5 inches. 

T his unfortunate detail of the Genah code is an historical one that is 
not likely to be altered in the future. 
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Summary: ichng Do these values need to be altered? 

Response: 0 The reported parameters are correct: 

Proceed with reconstruction. 

1 The reported parameters need to be altered. 

For any non-zero response to the parameter change query, the following 
set of requests are made: 

Summary: alat Enter the axial lattice spacing, "alat": 

Response: “a” Axial lattice spacing in units “units” 

Summary: drhx Enter the standoff distance, "drhx": 

Response: “b” The distance, in units “units”, from the hologram 

surface to the desired reconstruction surface. 

Summary: radius Enter the reconstruction radius: 

Response: “c” The radius of the desired reconstruction surface, 

in units “units” 

5.3 Output file name 

When an algorithm is constantly repeated or reused, the issue of filenaming is 
never a simple one. Genah gives almost complete flexibility in naming flies 
with two caveats. First, the name, including the full path, must be less than 
80 characters (see Section 5.1.2 for examples in Unix and VMS). Second, 
Genah does try to control the three-character extension to the fllename. 
This is most effective in normal reconstructions (Processing Option 0), where 
four output files are created. The same path and flle name are used, but the 
extensions begin with different letters depending on the data contained in the 
flle. These letters and their corresponding flle types are listed below in Table 
5.1. The last two letters of the extension are either two numerals or the letters 
“vw”. If a flxed wavenumber Alter is used for the entire reconstruction set, 
the two numerals are the cutofi" wavenumber multiplied by 10. For example, if 
kr(cutoff) = 6-0, the extension numerals would be “60”. If, on the other hand, 
a frequency-dependent Alter flle was used, the letters “vw” are appended as 
an acronym for “variable tyindow”. 
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Extension 

Data Type 

0 

ASCII processing description output file 

P 

Spatial domain Pressure, P{z, 0; r') 

V 

Spatial domain Velocity, V{z, 4>\ r') 

Order-limited Velocity, V( 2 ,^'; r') {n = ..., Umax) 

I 

Spatial domain Intensity, 7(z, r’) 

K 

Wavenumber domain Pressure without propagation, p(A:r, n; vq) 

G 

Wavenumber domain Velocity with propagation, v(A;^, n; r') 

N 

Axially reconstructed Velocity, V{z, n\ r') 

Q 

Axially transformed Pressure, P{k^, <f)-, tq) 


Table 5.1: First letter of filename extensions for Genah output files. 


Summary: fnamel Enter the Output File Name: 

Response: name.*mn mn = 10m.n = ^^K(cutoff) 

name. *vw “vw” for variable window (frequency-dependent 
wavenumber filter) 

(Note that the character is literal here) 
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Chapter 6 

Input and Output Files 


Perhaps the most cumbersome hurdles in the use of Genah are the format 
and structure of the data files for input and output. Again, it is the historical 
evolution of the code which led to the current situation. 

When it was originally developed on a PDP-11 platform for planar holog¬ 
raphy, the available storage (the 8-inch floppy) was quite dear, and the data 
needed to be stored as compactly as possible. This required the use of binary 
storage (4 bytes per real number, 8 bytes per complex number). When the 
broadband system was introduced, analyses required being able to access any 
particular hologram immediately, without having to step through the entire 
file each time. This required the direct access mode of storage. Direct access 
modes in FORTRAN-77 require a fixed record length. This immediately led 
to a conflict, in that different size holograms would demand different record 
lengths, and not all computer systems store record length information exter¬ 
nally. Since almost all holograms at that time had meshes of 64 by 64 points, 
it was decided that each scan line should constitute a record. 64 complex 
values require 512 bytes of storage, and this was chosen as the record length 
for all subsequent holograms. 

Thus, every hologram taken since 1985 at NRL is stored in direct access 
mode with a fixed record length of 128 longwords (4 bytes per longword). 
While complicating matters for different hologram sizes, this has made all 
previous experiments accessible to current software without modification. 
Also, both input and output files for every stage of processing share this 
same structure and format, so the same I/O software is reusable at every 
stage of processing. The format of these files is defined formally in Appendix 


45 




A. It is discussed here briefly so that the user can proceed with reading and 
writing holograms. 

A hologram is composed of a set of scan lines. For scan lines with M = 
2,4,8,16 or 32 axial locations, N = 64/M scans fit in each record. For M = 
64 locations, each scan line constitutes a single data record. For scans longer 
than 64 locations, additional records are used for each scan line. So, for 
instance, a 128-point scan requires two records; a 256-point scan requires 
four records. In short, the number of records required for each scan is the 
number of points in the scan divided by 64, rounded up to the nearest integer. 
The current version of Genah is limited to 256 axial locations per scan line. 

An individual hologram (at one frequency) requires as many records as 
there are scan lines times the number of records per scan line. For example, 
a typical cylindrical hologram of 128 points per scan line and 64 scan lines 
would require 128 records of 512 bytes each, for a total of 64 kilobytes of 
storage. A hologram with 32 points per scan line and 32 scan lines would 
require 16 records of 512 bytes each, for a total of 8 kilobytes of storage. 

The flrst record (512 bytes) of each file (whether input or output) is 
dedicated to a header record which contains “meta-data”, i.e., data about the 
holograms. This header is difficult to extract (again for historical reasons). 
It is critically important however, because it contains information on the size 
of the holograms, their number, and all of the pertinent information about 
the experiment needed by Genah for processing. The structure and contents 
of the header record are discussed in complete detail in Appendix B. Since 
this header record is always the first record, accessing a particular hologram 
record requires that 1 be added to the computed count to point to a given 
record in a specific hologram. 

It is implied above that a scan line may have any number of points. 
However, as currently written, Genah strongly restricts the allowable record 
sizes for reading and writing holograms and reconstructions. This originally 
grew from the fact that the number of points needed to be a power of two 
for the Radix-2 FFT routine used in Genah . Currently, Genah will only 
allow input files to have 16, 32, 64, 128 or 256 points per scan line exactly. 
If the data contain some other number of points per scan line, they must be 
preprocessed and rewritten, padded out to the next highest of these values. 
This rather restrictive practice is currently under review, and a future release 
of Genah will most likely provide much greater flexibility in the number of 
points per scan line. Even then, however, the data will be padded out to the 
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next power of two for the FFT. E the user needs even more flexibility, some 
other FFT routine will need to be substituted for the one provided in this 
Genah release. In any event, with the source code in hand, the user is free 
to alter the I/O to suit local needs. 

The available output sizes from Genah are also restricted. (This, too, 
is under review for a future release.) For spatial reconstructions (Processing 
Options 0, 3 and 5, see Section 4), only 16, 32, 64 and 128-point outputs are 
currently allowed. If the axial data have been zero-padded (izadd set to 1), 
then only the center M/2 points of an M-point scan line are written, so that 
the output line length is equal to the input line length. For wavenumber- 
domain outputs (Processing Options 1, 2 and 4), up to 256 points can be 
written out. In this case, it does not matter if the data had been zero-padded 
(up to a maximum of 256 points per line); all of the wavenumbers are written. 
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Part II 

Reference Manual 
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Chapter 7 

Theory of Nearfield Acoustical 
Holography 


7.1 Introduction 

The reconstruction of wave fields from holographic measurements has a long 
history in this century. It derives fundamentally from applying Green’s The¬ 
orem to the Helmholtz Equation for linear wave fields. While the optical 
version of holography is by far more well known, applications to acoustical 
fields have been known for decades. 

The idea that an acoustical field could be reconstructed at its source with 
a resolution that exceeds the acoustical wavelength limit dates to the seminal 
paper of Williams and Maynard in 1980 [1]. It was shown there that under 
certain circumstances, non-propagating (or evanescent) waves whose wave¬ 
lengths were smaller than the acoustic wavelength at a given frequency could 
be reconstructed. This reveals wave patterns distributed over the source sur¬ 
face which do not necessarily propagate to the acoustic farfield but which 
are critical to an understanding of wave motion on and energy propagation 
through the source structure. Thus, Nearfield Acoustical Holography (NAH) 
was born. 

What follows is a brief overview of the theory of NAH. For more detail, 
the reader is referred to the very thorough development in Maynard, Williams 
and Lee [2]. This work is quite complete and fieshes out many of the details 
glossed over here. The extension of the theory in Cartesian coordinates to 
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cylindrical coordinates is touched upon in that paper, but it is more elabo¬ 
rately treated in Williams, Dardy and Washburn [3]. The development of an 
algorithm for computer implementation of NAH is straightforward, but the 
treatment of errors and other issues is thoroughly covered in Veronesi and 
Maynard [4], a companion article to [2]. While the literature is rich with 
other articles (mostly dealing with implementations of NAH), these three 
form a useful core for understanding nearfield holography, and all three are 
found in the Journal of the Acoustical Society of America, which is readily 
available in most technical libraries. 


7.2 Analytical Development of Holography 

7.2.1 An Informal Sketch of Holography 

We begin our discussion of NAH theory with a simple motivation that reveals 
the concept behind reconstruction. This is followed by a more detailed anal¬ 
ysis of the functions involved for Cartesian Coordinates. Then, in Section 
7.4, the analysis is extended to cylindrical coordinates. 

Imagine a Cartesian half-space, that is, a three dimensional space extend¬ 
ing upward above a planar “floor” of infinite extent. Locations on the floor 
are denoted by coordinates x and y, and distance above the plane is mea¬ 
sured along coordinate z. Now imagine that this half-space is filled with a 
fluid which will bear acoustic disturbances, but that no source of disturbance 
is located anywhere above the floor. The propagation of acoustic pressure 
waves in the half-space is described by the wave equation: 

1 

V^p(a;, y, z; t) - y, z; t) = 0 , z > 0 , (7.1) 

where c is the speed of sound in the fluid. Assume the time variation of 
any waves in the space (or on the surface z = 0) is harmonic, that is, of the 
form with some frequency oj. If we let p(x, y, z; t) = p(x, y, z)e“’“'* 

and substitute this into 7.1, we arrive at the Helmholtz equation: 

V^p(x,i/,z)+k^p(x,p,z) = 0 , z>0,P=(jy. (7.2) 

In essence, we have performed a Fourier transform on the wave equation to 
reduce the dependence on time to a set of independent equations for difiFerent 
(arbitrary) frequencies. 
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We formalize the concept of the Fourier transform by defining the trans¬ 
form pair for a spatial variable 

/ OO 

pK)e-'W^, (7.3) 

*oo 

An important property of this transform is that when it operates on the 
second derivative of a function, the result is just the negative of the squared 
frequency times the transform of the function itself: 

/ OO 

=-klP(k). 

The operator in Cartesian coordinates is simply the sum of the second 
derivatives in each direction. 

If we apply the Fourier transform to the x and y directions of 7.2, this 
property of Fourier transforms yields 

S 

-klP(k„ky,z)-klP(k,,ky,z)-¥-^P(k,,k,,z)+k‘‘P{k,,k,,z) = 0 ,z > 0 , 

(7.5) 

where k = uj/c is known as the propagating wavenumber in the fluid. It is 
related to the fluid wavelength A at frequency u by the equation k = 27r/A. 
When the terms in equation 7.5 are regrouped, we obtain 

cP 

ky, z) + [k^ - {kl -I- kl)]P{kx, ky,z) = 0 , z>0 . (7.6) 

This expression is a simple, second order, ordinary differential equation whose 
solution is well known; 

P{k^, ky, z) = . (7.7) 

Note that by expressing certain quantities in their related Fourier domains 
through repeated applications of the Fourier transform, we have moved from 
the wave equation, which is very difficult to solve, to a simple ODE. 

Because of the way we have defined the coordinate system and the re¬ 
lation between time and space in the harmonic functions, the first term in 
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7.7 represents outward travelling waves (up from the floor) and the second 
term represents inward travelling waves (down toward the floor). However, 
the statement of the problem deflned the half-space above the floor as be¬ 
ing source free; no waves could be coming downward from above the floor. 
This boundary condition forces 5 = 0 in 7.7. Application of this source free 
condition is crucial to the success of NAH. Without it, we can’t tell unam¬ 
biguously whether a measured acoustic disturbance is an outgoing wave or 
an incoming wave. 

To solve our equation, we need another boundary condition. The floor, 
or source surface, provides that condition. The pattern of waves on that 
surface is the only available source of acoustic energy in the fluid. The 
acoustic pressure right on that surface has some distribution over x and y 
that we describe via the spatial Fourier transform on those coordinates, that 

oo 

p(x. S'; ^ = 0) = ^ / /~ ^ dk, . 

But P{kx,ky]z = 0) is precisely the coefficient A of equation 7.7 (set z =■ 0 
in this equation to see this). Thus, defining Po{kx, ky) = P{kx, ky\z = 0), we 
obtain an expression for the wavenumber domain pressure everywhere above 
the floor: 

P{kx,ky;z) = Po{kx,ky)e''^^^ (7.8) 

Applying the inverse Fourier transform again, we obtain 

pix,y,z) = ^l j"^Po{kx,ky)e''^^^P^^ye'^^^dkx dky . 

This equation immediately implies that the distribution of waves on the sur¬ 
face bounding the half-space completely determines the pressure throughout 
the half-space. This is the key to holography, which literally translates to 
“whole writing”; write or graph the solution on the boundary, and the whole 
field is determinable. 

Through all of this explication, we haven’t discussed the nature of the 
wavenumber k^. If a plane, harmonic pressure wave is propagating into the 
half-space in some direction away from the floor, we describe its direction 
with the vector k = (fcjej,,A;^e^) where (ex,ey,e^) are the unit vectors 
in x,y and z. The size of k is determined in the usual Euclidean sense, 
k"^ = k‘^ + ky + k"^. Since k is known for a given fluid and frequency, and kx 
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and ky are known from the distribution on the bounding surface, it is easy 
to see that is definable in terms of them: 

k't = k^- {kl + kl) 

Although k is fixed by the fluid and frequency, no such restrictions are placed 
on kx or ky\ we can impose (or imagine) any wave distribution we desire on 
the bounding surface. This means that the quantity can be negative as 
well as positive, and that its root k^ can be either real or imaginary. Herein 
lies the key to nearfield holography. 

Suppose that a harmonic wave of frequency u exists on the floor with a 
wavelength much shorter than that supportable in the fluid at that frequency. 
It propagates across the floor with a trace wavenumber k^ = k^ ky. This 
trace wavenumber is greater than the fluid wavenumber k, so k^ becomes 
imaginary. Then the solution to equation 7.8 becomes 

P{k„k„z) = Po{k„k,)e‘^‘'^-> 

= P^(K,k,)e-'‘-. 

This evanescent wave distribution on the floor decays exponentially in the z 
direction. This is the definition of an acoustical nearfield: wave distributions 
on a source cannot propagate to the farfield. Solutions to equation 7.8 must 
then be classified into one of two categories depending on the relationship of 
the trace wavenumber on the bounding surface and the fluid wavenumber at 
the same frequency: 

P{kx,ky,z) = PQ{kx,ky)e'^‘’' , kt<k\ (7.9) 

= Po{kx,ky)e~^‘^ , kt> k . (7.10) 

The power of nearfield holography lies in the ability to describe the entire 
wave distribution on the source surface, including both evanescent and prop¬ 
agating contributions. 

So how does all of this apply to measuring a hologram and making recon¬ 
structions? In equation 7.9, suppose the variable z was replaced by a specific 
value Zff, the distance to a hologram surface. Further, imagine that the pres¬ 
sure distribution was measured over all x and y at this surface, then Fourier 
transformed. By inverting equation 7.9 directly, the pressure distribution at 
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the source surface 2 = 0 would be entirely characterized: 

Po{k,,ky) = kt<k (7.11) 

= P{k^,ky,ZH)e^‘^^ kt>k. (7.12) 

This is the informal definition of a reconstruction by back propagation. Note, 
however, that inversion is not a free lunch. Any such inverted system is an ill- 
posed problem, and questions of uniqueness of the solution are not trivial, nor 
can they be ignored. While we only touch on these issues in this document, 
it is easy enough to see what happens to Pq in the case kt > k when zjf is 
large in relation to the wavelength A^. Since, in experiments, most of the 
measurement noise is manifested in the highest wavenumbers (where kt k), 
reconstructions can quickly become meaningless if 2 ^/ is chosen too large and 
appropriate wavenumber filtering is not applied. 

7.2.2 A Formal Derivation of Holography 

The following derivation is lifted mostly from [2]; the reader is encouraged 
to study that reference more closely. 

We first give some definitions. We assume that some source is creating 
an acoustic pressure field p(r, t) in a space spanned by r over a time t which 
satisfies the homogeneous wave equation 7.1. We assume that some surface 
S bounds this space, and that there exists a Green’s function G(r|r 5 ) which 
satisfies the homogeneous Helmholtz equation 7.2 for r inside S and vanishing 
(or having a vanishing normal derivative) for r = r 5 on S. That part of 
S which is not at infinity is in practice a level surface of some separable 
coordinate system which is in close contact with the sources. In the sketch of 
Section 7.2.1 above, S was the “floor” of the semi-infinite half-space. Finally, 
we assume there is another surface H which is a level surface parallel to S 
for which either or its normal derivative can be measured. This is 

the hologram surface. 

First, the data collected as p{Tff,t) are transformed to the frequency 
domain: _ 

poo 

p{th,u})= p{TH,t)e^^ dt, (7.13) 

J — OO 

at which point the wave equation becomes the Helmholtz equation as shown 
previously. For a fluid in the space which will support acoustical waves with 
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a compressional speed c, a wave propagating with frequency oj is defined as 
having a wavenumber k = ojIc. In using the Helmholtz equation, we can deal 
with independent frequencies u. We are now concerned with solving for ^(r^) 
from p{Tff) at a given frequency u. Note that p{Tff) is a complex quantity, 
possessing both amplitude and phase (relative to a reference) information. 

We begin the spatial processing by examining the Helmholtz Integral 
Equation. This equation derives from application of Green’s Theorem to the 
Helmholtz equation, in this case for acoustic pressure. In its most general 
form, it has two terms: 

“ Mrs)|^(r|r5)) dS (7.14) 

where G(r|r 5 ) is the Green’s function satisfying the Helmholtz equation 
which propagates waves from ts to r in the space bounded by S, and d/dn 
is the derivative in the direction normal to dS. For holography, we construct 
G(r|r 5 ) as the sum of the free-space Green’s function plus another term 
which will force G to zero on the boundary ts. This will then satisfy the 
homogeneous Dirichlet conditions on 5, that is, G(r|r 5 ) = 0. This forces the 
first term in 7.14 to vanish, leaving only 

Note that the integration is over the surface S, which we take to be the first 
two spatial coordinates described by r^. If 5 is a level surface of a separable 
coordinate system, and we fix that level by taking the third coordinate to be 
constant, then G can be found analytically and equation 7.15 becomes a two 
dimensional convolution integral. 

In holography, we do not know the field on the surface S. We know it on 
another level surface H. Using generalized coordinates ^i, ^ 2 , and ^ 3 , we can 
express 7.15 as 

-e;,e 2 > (7.16) 

where the primed coordinates vary over S and the unprimed coordinates vary 
over H. In this expression, C = is a constant (we know the distance 

from our holography surface to our source surface). In principle, this is a 
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two dimensional convolution, as mentioned above. However, we know the 
pressure over H\ how can we extract the pressure over S from inside the 
convolution? 

The equation is invertible through the use of the convolution theorem. 
Without going into detail, the convolution theorem states that the product 
of the Fourier transforms of two functions is equal to the Fourier transform 
of the convolution of the functions themselves, that is 

/ OO 

^{x)4>{x - y)dy = > 

•OO 

where T is the Fourier transform operator, and = T‘ip{x) and = 
If we take the two dimensional Fourier transform in and ^2 of 
equation 7.16, we get 

Pm = P{^i)G'{0 (7.17) 

where G' is the Fourier transform of the Dirichlet Green’s function relating 
the transformed surfaces by their separation distance, (. The other vari¬ 
ables have been suppressed for clarity. We can now (informally) invert this 
expression and solve for the unknown P{^§) in terms of the known Pi^^): 

P((l) = lP((f)G'-'(0] (7.18) 

This equation directly relates the source distribution to the distribution 

measured on H, but we don’t yet know how to propagate the measurement 
to other surfaces ^3 = constant. To do this, we rewrite equation 7.17 for the 
new surface ^ 3 : 

P(ii) = P{4)&(U-il)- (7.19) 

We have an expression for P{^i) from equation 7.18, and we can substitute 
it into 7.19. Doing so, and taking the inverse Fourier transform over the first 
two coordinates, we arrive at the general expression for holographic recon¬ 
structions throughout the space bounded by S in terms of the measurement 
taken at H: 

= . ( 7 . 20 ) 

We have arrived at the fundamental expression for nearfield holography. 
It tells us what we need to do (theoretically, at least) to determine the dis¬ 
tribution of waves on some separable surface. First, measure the waves on a 
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surface parallel to the source surface (and close to it). Fourier transform the 
measured distribution over the surface. Develop the Fourier transform for the 
Dirichlet Green’s function separating the hologram and source surfaces and 
the function separating the reconstruction and source surfaces. Invert the 
former and multiply both with the measured function for all wavenumbers. 
Inverse Fourier transform the result. Repeat this process at each individual 
frequency of interest. A set of reconstructions ensues. 

We stress here that this inversion is informal. There is an entire body of 
mathematics literature which holds that formation of is at best risky, 
at worst impossible, and that the entire enterprise of reconstruction does 
not form a unique solution. We do not intend to delve into the details of 
“the Inverse Problem” here. Suffice it to say that the inversion works in an 
engineering sense, in that under the controlled circumstances of real NAH 
experiments and processing, reconstructions are similar to the “true” answer 
to within some error bounds which are tolerably small. 


7.3 Holography in Cartesian Coordinates 

We have, in essence, solved the problem of nearfield holography analytically. 
It remains to derive the Green’s functions and their Fourier transforms for 
the coordinate systems of interest. 

In Cartesian coordinates, the separable surfaces are planes. As in the 
informal sketch presented in Section 7.2.1, we align our axes so that the 
surfaces are planes in (x, y), separated by some distance zu — zs- The Green’s 
function which satisfies the homogeneous Dirichlet boundary condition on zs 
is given by 


G{x,y,z\x' ,y' ,z') = 

exp [i kyj{x — x'Y + {y — l/Y + (-2^ “ ^'Y\ 
ylix- x'Y + {y- y'y + {z- z'f 
exp[ik^J{x - x'Y + (y - iff + {z- z' - 2zsf] 
y'(x - x'Y + (t/ - ?/')2 + (z - - 2zsf- 

(We do not develop this function here: it can be found in Morse and Feshbach 
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[5].) Fixing z' = zs and taking the normal derivative we obtain 

( \ t ' \ 

—{x,y,z\x,y,zs) = 


^ d , exp[ik^{x - x'f + {y- y'f + v, 

5 a V _ x'Y + {y- iff + a2 (--s) ‘ 


These expressions apply, of course, to the spatial domain. To apply them 
in our holography algorithm, we require the Fourier transform of 5G/5n. 
Fortunately, it can be found analytically: 

G'ih,ky,z) = kl + kl<k\ 

= , kl + kl> k^ . (7.23) 

The wavenumber quantities in the exponentials are the effective z wavenum¬ 
bers. As in our informal sketch of holography, the issue of the trace wave- 
number enters and gives rise to two exclusive solutions. As long as the trace 
wavenumber k^ = kl + ky is less than the fluid wavenumber k^ = {ufcy at 
frequency u, the waves distributed on surface S will propagate into the fluid. 
When the trace wavenumber exceeds the fluid wavenumber (the wavelength 
on the surface is less than that in the fluid), the wave on S decays exponen¬ 
tially along 2 : into the fluid. In a two dimensional image of the wavenumber 
domain, the boundary formed by this relationship is the “radiation circle”; 
its radius is the trace wavenumber, kf. This radius increases linearly with 
frequency (since the fluid wavenumber is proportional to frequency), forming 
a “cone” in wavenumber-frequency space. 

Now that we are in possession of the proper expression for the Dirichlet 
Green’s function, we can express the holographic reconstruction relation for 
planar coordinates. Substituting equation 7.23 into equation 7.20 yields 

p{x,y,z) = , 

= kl + kl> .(7M) 

Note that zs does not explicitly appear in this expression. This is due to the 
simplicity of the product &~^(zh — zs)G'{z — zs) in Cartesian coordinates; 
the term with zs factors out. More generally, the surface S only establishes 
the range of validity of the final expression. While it could certainly be the 
z to which we reconstruct the data (and usually is), as the source plane, it 
does not play an explicit role in propagating waves from zh to z. 
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7.4 Holography in Cylindrical Coordinates 

Cylindrical coordinates are described in terms of axial distance, a circum¬ 
ferential angle 0 about the origin axis and a radius from that axis of r. A 
bounding surface S consists of an infinitely long, right circular cylinder cen¬ 
tered on the 2 :-axis of radius rs', only the space r > rg is describable by 
holographic reconstructions. A hologram measurement surface H would be 
over another right circular cylinder concentric to S at radius rg. 

In this coordinate system, the wavenumber domain co nsists of the ax¬ 
ial wavenumber, n, the circumferential order, and , the radial 

wavenumber. Note that kr is analogous to the k^ of Cartesian coordinates 
in that it is derived from the fluid wavenumber and a trace wavenumber. 
This implies that it, too, can be either real and propagating (when k^ < k) 
or imaginary and evanescent (when k^ > k), as we will see momentarily in 
the solution. Note, too, that kr does not depend on the circumferential or¬ 
der n. This does not imply that circumferentially oriented waves cannot be 
evanescent Zi/ie; we will be able to say more about this when we examine the 
solution. Finally, note that while k^ and k are continuous variables (and kr 
by extension), n is drawn only from the set of integers. To be continuous, 
a wave solution around the circumference of a circle must have an integer 
number of wavelengths. (Zero is a possible solution.) This has very strong 
implications for the use of discrete transforms in the algorithm for nearfield 
holography, as will be discussed in Chapter 8. 

We take a slightly different route in developing the expressions for cylin¬ 
drical coordinates than that employed for Cartesian coordinates. (This de¬ 
velopment is lifted from both [2] and [3].) Fixing the radius at r^, we expand 
the pressure distribution there in terms of (orthogonal) cylindrical wavefunc- 
tions: 


1 °° roo 

The wavefunction coefiicients Pn{kz\rtj) are found by taking the two di¬ 
mensional Fourier transform of equation 7.25 and invoking the orthogonality 
conditions: 

1 /•2)r r<x> 

Pn{kz\r]j) = ^ j dk^)d^ . (7.26) 
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Because p{z, (j)\ r) was expanded in orthogonal wavefunctions, equations 7.25 
and 7.26 can be seen to be a Fourier transform pair (noting that the inverse 
transform on n, which is discrete, becomes a Fourier series). In other words, 

p{z,4>\rH) = J^~^[Pn{kz\rn)] and 

Pn{K\rH) = T\p{zA\r„)] . (7.27) 

We now have expressions relating pressure wave distributions to the coeffi¬ 
cients of cylindrical wavefunctions via the Fourier transform. 

We now need to find how to propagate a known distribution radially 
from one surface to another. In equation 7.25, we fixed and expanded 
in wavefunctions only over a cylinder of radius r^. The general solution in 
cylindrical coordinates contains radial functions as well which are expressed 
in terms of the Hankel functions, = J„(x) ± iYn{x), where Jnix) 

and y„(x) are the Bessel functions of order n and argument x of the first and 
second kinds, respectively. Since the holographic space is entirely exterior to 
the bounding surface S and is source free, only outgoing waves are allowed. 
This restricts the solution to H^\x), the Hankel function of the first kind 
of order n and argument x, which we will refer to here as simply Hn{x). 
Equation 7.25 generalizes to 

p{z,<f>,r) = — (/ Fn(W„(fc.r)e‘^'^dA:,)e'”^ (7.28) 


where, as mentioned above, When k^ > k, the argument to the 

Hankel function, /c^r, becomes imaginary, and the Hankel function becomes 
a Modified Hankel function KniKr). Whereas the Hankel function behaves 
asymptotically as exp{ikrr) (with cylindrical spreading), the Modified Hankel 
function behaves as exp{—krr); when the radial wavenumber becomes imag¬ 
inary, the axial trace waves no longer propagate, but become evanescent, as 
in the Cartesian case. 

The radial wavefunction in equation 7.28 forms the Green’s function for 
this coordinate system. It does not participate in the Fourier transforms on 
z and (j). We can be rewrite 7.28 for any radius, including rg: 


1 fOO 

p{z, (t>.r„) = — Y. i Pn{kz) Hr^iKTH) dk,) . 
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Note that the P„ are not functions of r, but eigenvalues of the cylindrical 
wavefunction expansion. The radial (propagation) information is contained 
in the Green’s function, HniKr). In other words, the pair of relations 7.27 
held only for a fixed value of r. To relate distributions at different radii, 
we must express 7.27 with the radial propagator included. Doing this at 
an arbitrary reconstruction radius r and at the specific hologram surface rn 
yields 


p{z,(l),r) = T ^[Hnikrr)Pn{h)] and 
p{z,(l),rn) = P'~^[HnikrrH)Pn{kz)] ■ 

Solving for the coefficients Pnikz) in the second expression and substi¬ 
tuting into the first expression, we arrive at the relation for holography in 
cylindrical coordinates: 

p(z,4>,t) = T~'[mKr)H~'(k,rH)P„{k,)] \ k, < k , 

= T-'lK„(Kr)K;;'{krrK)P„(k,)]; k, > k . (7.29) 

These equations are valid for any r > rs, the smallest cylinder which will 
enclose all the sources. Compare this expression with that of equation 7.24 
in the previous section: only the Green’s function has changed. 

As an aside, we consider the evanescenthTie behavior of circumferential 
waves smaller than the fluid wavelength at a given frequency, as mentioned 
above. The higher the order n for a circumferential wave, the smaller the 
wavelength on a given cylinder, since = 27rr/n. When n becomes large 
relative to the argument k^r, the Hankel function can be asymptotically 
approximated by 

When xjn < 1, the second term dominates, and the Hankel function decays 
as the nth power of radius. At a radius tq = An/27r, the first term begins 
to dominate, and the power law decay gives way to the cylindrical spreading 
farfield proportional to l/-\/r. This is demonstrated clearly in Section I.A.l 
of [3]. 
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7.5 Other Field Quantities 

So far, we have concentrated on reconstruction using the Dirichlet Green’s 
functions. In acoustics, this translates to the pressure field. However, the 
power of nearfield holography is greatly extended by the introduction of re¬ 
constructed velocity. In fact, normal velocity on the surfaces of cylindrical 
shells is probably the most important working quantity for structural acous¬ 
tics research. 

In addition, pressure and velocity reconstructions can be combined to pro¬ 
vide acoustic intensity. This quantity is the key to source and sink location on 
radiating structures as well as computation of power flow and total radiated 
power. In many cases, this last quantity is used as a metric for comparing dif¬ 
ferent cases, especially simulations using coupled Finite/Infinite/Boundary 
Element codes. 

We won’t provide detailed derivations of reconstructed velocity; instead, 
we simply state the formula needed to compute them. 

7.5.1 Reconstructed Normal Velocity 

There are two paths which will lead us to a reconstruction formula for nor¬ 
mal velocity: through Newton’s law and through the Neumann boundary 
conditions. They are, of course, different views of the same problem. We 
approach the topic through Newton’s second law. 

In analogy with the statement F = ma, we relate the product of the time 
rate of change of the fluid velocity (acceleration) and the density (mass) to 
the gradient of the pressure (force): 

(7.31) 

where each component of dv/dt is related to the affiliated component of the 
gradient operator. Assuming harmonic time dependence, we can directly 
relate the velocity to the pressure: 

V = T—Vp . 
lujp 

This simple expression gives us a formula for deriving the operator for trans¬ 
forming from measured pressure to reconstructed velocity. Since we will be 


(7.32) 
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operating in the wavenumber domain, we Fourier transform 7.32 and substi¬ 
tute the propagated pressure as shown in 7.19. This yields 

V,(6,6.&) = , (7.33) 

LUJ ij 

where rj represents the individual vector components. For radiation, the 
component in the direction normal to the source surface (the propagation 
direction) is of greatest interest. Since P is not a function of the normal 
coordinate, the normal gradient does not operate on it, but only on the 
Green’s function propagator. Of course, the work here involves computing 
the transform of the gradient operator, then applying it to the propagator. 
More details may be found in [2] and [3]; we simply state the results. 

In Cartesian coordinates, the normal velocity propagator appears only 
slightly modified from the pressure propagator of equation 7.23: 

G',{K, ky,z - zh ) = (y) ; kl + k^K , 

= (^) ; kl + kl> k\ (7.34) 

In cylindrical coordinates, the normal component is in the radial direction. 
In this case, the propagator of equation 7.29 becomes 

G'^{kz,n,r- th) = -K [JT'(fc,.r)i7"^(A:,.7'ff)] ; k^ < k , 

= -kr [K{Kr)K;:\krrff)] ; k^ > k , (7.35) 

where the prime on the numerators indicates the derivative of the function 
with respect to argument. Since the derivatives of Bessel functions are eas¬ 
ily determined from the Bessel functions themselves via a recursion relation, 
these propagators are no more difficult to compute than the pressure propa¬ 
gators. 

7.5.2 Acoustic Intensity and Radiated Power 

With the pressure and normal velocity distributions available over a surface, 
the power flow per unit area, or intensity, is computable. The fundamental 
formula is simply 

I = iRe(pV) , (7.36) 
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where v* denotes the complex conjugate of the velocity and Re represents 
the real part of the complex product. This quantity is particularly valuable 
in the nearfield, where it can map power flow. By definition, in the farfield 
all power flow is outward. 

Once the intensity is computed everywhere over the bounding surface, 
the total radiated power can be computed for each frequency. By definition, 
the power is the area integral of the intensity: 

n = i / Re{p\*) dS . 

^ J s 

For the planar system, dS = dxdy\ for a cylindrical system, dS = rd(j)dz. 
Note that Parseval’s theorem guarantees that we could also compute total 
power from the integral over wavenumber of the modulus squared pressure 
coefficients in the wavenumber domain. Such an integral would be limited to 
the supersonic wavenumbers, since they are the only ones that radiate. 


(7.37) 
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Chapter 8 

Implementation of the NAH 
Algorithm 


In Section 7.2.2, we developed the basic recipe for reconstructing nearfield 
acoustical holograms. We transform pressure time series taken at a mesh of 
locations describing a separable coordinate surface to the frequency domain, 
then rearrange the data into holograms at each frequency. To perform recon¬ 
structions, we Fourier transform the measured distribution over that surface 
to the wavenumber domain. With the Fourier transforms of the Dirichlet 
Green’s functions for the hologram and reconstruction surfaces in hand, we 
multiply the wavenumber domain hologram data by the Green’s functions 
point for point over all the wavenumbers. We then inverse Fourier trans¬ 
form the result back to the spatial domain. We repeat this process at each 
individual frequency of interest and develop a set of reconstructions. 

To accomplish this set of operations on a computer, every aspect of the 
process must be discretized. Discretization is straightforward, but it presents 
a set of challenges about which the experienced holographer must be aware. 
The side effects of discretizing operations which are modelled continuously 
impact not only the processing of holographic data, but the design of experi¬ 
ments intended to collect the holograms in the first place. The reader is once 
again referred to [2] and [3], but the definitive discussion of algorithms and 
implementation is presented by Veronesi and Maynard [4]. 
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8.1 Transition to a Finite, Discrete Domain: 
Assumptions 

We begin with a discussion of the issues involved in discretization, followed 
by the actual formulae involved. We then examine how discretization impacts 
the experiments and the data in nearfield holography. 

The two most important issues in implementation of holography are the 
truncation of the infinite domain to a finite domain and the discretization 
of the continuum to some finite mesh. The development of the expressions 
for nearfield holography depended on the bounding surfaces being of infinite 
extent. In planar holography, this is the entire plane at some zs', in cylin¬ 
drical holography, this is the infinite cylinder at rs- In reality, we truncate 
this domain to some aperture of length L. In cylindrical holography, this 
forms a right circular cylinder of extent —L/2 < z < L/2. In order for the 
infinite integrals of the holography algorithm to remain accurate, we must 
assume that the function representing the measured quantity 'tp{z, 4>, r) van¬ 
ishes outside this domain. This might seem ridiculous in the free field outside 
a finite, radiating shell. However, two things come to our rescue. First, one 
usually takes holography data as close to the source as possible (more on this 
below). Since the acoustic pressure levels are typically quite high in the re¬ 
gion close to the radiator, they tend to fall below the measurement dynamic 
range within a relatively short distance (axially) from the body. Second, 
one usually “windows” the data (axially), forcing a smooth taper to zero 
value at the ends. Together, these conditions suffice to meet the vanishing 
function criterion. To help enforce the first condition, the aperture L must 
significantly exceed the length of the body under study. At NRL, we usually 
require a 50% overscan (i.e., 1/4L past each end of the shell) for cylindrical 
experiments. Aperture length will be discussed below in more detail. 

The discretization of the measured pressure field depends on two more 
assumptions in the spatial domain and a concomittant assumption in the 
wavenumber domain. First, we must assume that the continuous wave field 
can be replaced by a discrete, finite set of values which sufficiently represent 
the field. These values are assumed to be the average pressure over a patch 
spanning the discretization interval. Further, it is assumed that the continu¬ 
ous wavefield can be constructed from a set of point measurements, with one 
point per patch. This assumption strongly impacts the sampling interval] 


68 




the wave field cannot vary in a manner which is not captured by the discrete 
samples. This can be formalized with Shannon’s Sampling Theorem, whose 
primary result is commonly known as the Nyquist criterion: the sampling 
interval must be less than 1 /2 the wavelength of the shortest wave expected 
in the system. Otherwise, aliasing of higher wavenumbers into the sampled 
field will result. This leads to the third assumption, that a discrete, finite 
set of wavenumbers is sufficient to characterize the Fourier transform of the 
wavefield. This criterion impacts both the aperture length and the sample 
interval, and will be discussed below. 


8.2 Discretization of Integral Operations 

Based on the assumptions presented above, we are ready to examine the ef¬ 
fects of discretizing the continuous operations outlined in Chapter 7. The 
first of these is the two dimensional convolution known as the Rayleigh inte¬ 
gral, which allows the wavefield at one surface to be propagated to another 
surface using the appropriate Green’s function. We assume that the acoustic 
pressure field has been measured at a set of points describing a mesh within 
the measurement aperture, and that the size of this mesh is M x N. (In 
cylindrical holography, this would be M points along the axis and N points 
around the circumference.) Then the convolution integral becomes (see [4]) 

V-l M-l 

ip{p,q]r) = ^ p{m,n-,rH)G{p-m,q-n\r- th) . (8.1) 

n=0 m=0 

In this expression, p and q are the indices of the reconstructed point at radius 
r, m and n are the indices of the hologram data values at radius is the 
reconstructed quantity and G is the Green’s function propagator associated 
with Ip. Veronesi [4] goes to some length in describing methods of computing 
G; for our purposes, we simply sample it in the wavenumber domain as we 
now describe. 

The power (and speed) of Genah lies in performing the convolution of 
equation 8.1 in the Fourier domain, where it becomes an inner product. To 
do this, we use the Discrete Fourier Transform (DFT), implemented through 
the Fast Fourier Transform (FFT) algorithm. The DFT is formalized as 
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follows; 


1 N-l M-l 

Pn{kn.\r) = -^Y, e-''"'"""'' E ■/>,; r) AzA,/, , (8.2) 

9=0 p =0 
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Here, L is the length of the axial aperture, M is the number of points along 
that length, and N is the number of angles around the circumference of the 
measurement cylinder, as used above. Substituting for (l/27r)A^A(^ in 8.2, 
we obtain 


Pn{^m^ 


r N-l M-1 

^ •''' 9=0 p =0 


(8.3) 


We could conceivably compute the transform of G in the same way (again, 
refer to [4]). However, since analytical expressions for G are already available 
to us (see Sections 7.3 and 7.4), we find it much more convenient to simply 
sample those functions over the wavenumber set. Described by Veronesi [4] 
as Method 5, it is shown to be quite accurate for backpropagation over small 
distances. It is much less accurate in forward propagation over distances of 
more than a few wavelengths than some of the other methods he describes, 
but these instances are considered to be of far less interest to us. 


8.3 Mitigating Side Effects in Truncation and 
Discretization 

We now return to some of the issues brought up in examining the assumptions 
needed to truncate and discretize the hologram surface in Section 8.1. 
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First, we investigate the length of the measurement aperture, L. (In 
planar holography, both x and y require the same analysis.) One important 
criterion for L has already been described. This is the need for the wave field 
to vanish at the ends of the truncated aperture in order to satisfy the infinite 
convolution integrals. As mentioned, this requires that the measurement 
aperture overscan the shell by some amount (often about 25% past each 
end) to insure that the pressure falls to the bottom of the dynamic range. In 
addition (and often because this criterion fails, especially at low frequencies), 
we usually apply a tapered window to the data axially to ensure the vanishing 
constraint. Windowing as implemented in Genah is discussed in greater 
detail in Section 5.2.1. 

Another criterion in choosing the size of L is interference from replicated 
image sources. By their nature, DFT’s are periodic over a period 27rM/L. 
This implies that the Fourier transform representation of the hologram aper¬ 
ture repeats axially, and by analogy, the hologram itself repeats every L. If 
the hologram has significant contributions to the wave field near the ends of 
the aperture, the images of this source will contribute during the propagation 
process, interfering with the true propagation. Thus, L must be large enough 
to isolate the source and move its virtual images substantially away from the 
real aperture. Another way of looking at this is to consider the relationship 
between L and Ak^: 

Afc, = ^ . (8.4) 

If the source fills the aperture and L is not large enough, Ak^ will be too 
coarse to accurately represent the wave distributions; the field values at in¬ 
dividual km = mAk^ become contaminated by the interfering image sources 
and have the wrong amplitudes and phases. This effect can sometimes be 
mitigated by a combination of windowing and zero padding. In zero padding, 
the data are extended to the next power of 2 points with zeros, doubling the 
effective aperture and the distance to the nearest images. When windows are 
used to taper the wave field to zero at the edges to meet the padded aper¬ 
ture, effective control of truncation errors and replicated sources is gained. 
However, these signal processing “tricks” are not a panacea for good experi¬ 
mental design. The best way to guard against the effects of spilling over the 
aperture and having replicated sources is to take preliminary test scans of 
the source at the lowest and highest frequencies of interest and at frequencies 
where the source strength is large, and to analyze their impact at the ends 
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of the aperture. 

It is of greatest importance to note that these issues do not come into 
play for the circumferential direction. This dimension is, by definition, finite. 
The convolution integral on (j) is limited to the range (0,27r). Thus, there are 
no truncation or replicated image side effects in the (f> direction. 

Also of concern in designing a scan and later processing it is the sam¬ 
ple interval, Az. This interval directly determines the maximum available 
wavenumber in the axial direction, 

O'jr 

kJmax) = (8.5) 

As stated above, this interval must be sufficient to capture the smallest wave¬ 
lengths likely to occur in the system. (As an aside, these are almost always 
the axial flexural waves on a cylindrical shell at the highest frequency of in¬ 
terest.) This choice is bounded by the Nyquist interval, 1/2 the minimum 
wavelength. In practice, this upper bound should be given wide berth: ^Xmin 
to jAmin is not unreasonable. In our controlled laboratory setting, however, 
we have been able to approach very close to the Nyquist limit and still get 
decent reconstructions. 

Although one would like to go as high in wavenumber as possible, it must 
be remembered that the highest wavenumber contributions, almost always 
evanescent (or, equivalently, subsonic), grow the most rapidly in back propa¬ 
gation. This results in the destruction of the reconstruction’s dynamic range. 
The solution to this is judicious use of wavenumber filtering, as described in 
Section 5.2.4. Nonetheless, there is no need to push k^f^max) so high that the 
Alters cut off" well below it for all frequencies. This equivalently small 
would mean that the system was too oversampled, a waste of very expensive 
resources. The additional errors introduced at high k^ due to discretization 
are treated in great detail in [4]: they are not discussed here. 

The same principles of sampling do apply to the circumferential direction. 
Although the “aperture” (27r) is fixed, the number of angles sampled around 
the circumference must be sufficient to capture the smallest wavelengths 
expected in the system. 

Finally, one is not, unfortunately, free to choose L and Az independently 
of one another. Because of the Radix-2 FFT’s used in Genah , M and N 
must each be a power of 2. In the case of M = LfAz, the choice of either L 
or Az constrains the other quantity to within a power of 2. This often calls 
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for some compromise. Our priority is usually placed on Az, since many of 
the effects of the truncated aperture can be mitigated by windowing and zero 
padding, as described above. Most of our scans have been either 64 or 128 
points in length; 32 is seldom enough for meaningful shells and 256 is often 
too expensive in terms of acquisition time. We have almost always used 64 
angles circumferentially due to the use of point drives in the shells; in cases 
of ring drives or axisymmetric scattering, 32 or 16 angles sometimes suffice. 


8.4 Resolution, Dynamic Range and Numer¬ 
ical Precision 

We have mentioned several times previously the need for and the positive 
effects of locating the hologram measurement surface as close as possible 
to the source surface. This section will quantify that separation distance 
{d — zjj — zs in Cartesian coordinates, d — rg—rs in cylindrical coordinates). 
It turns out to be directly related to the achievable wavelength resolution 
and to the measurement dynamic range. Thus, the value of d plays a part 
in choosing Az and L as well. 

Assume that a holographic measurement system has a given dynamic 
range D, measured in decibels (i.e., 20log{Pr„axfPmin))- For a digital acqui¬ 
sition system, the theoretical best-case baseline is based on the number of 
bits in the A/D converter. Accounting for the fact that the acoustic pressure 
has both positive and negative values (thus giving up one bit to the sign), 
a 10-bit converter gives a precision of 1 part in 512 (about 54 dB), a 12-bit 
converter gives 1 part in 2048 (about 66 dB) and a 16-bit converter gives 1 
part in 32768 (about 90 dB). Of course, transducer, preamplifier, amplifier 
and timing/positioning/phase noise all reduce this theoretical value, but they 
differ for any two installations. Hopefully, their contributions can be kept to 
less than 12 dB (the two least significant bits). 

For a given evanescent wave component at wavenumber A;(., the field will 
decay exponentially as the transducer is moved away from the source surface. 
At some distance, the component’s amplitude will fall below the minimum 
threshold of the measurement system. This is how dynamic range relates to 
physical distance: 

jq-O/20 ^ 
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The achievable spatial resolution is bounded by the Nyquist rule, as dis¬ 
cussed above: 

r, 27r 

Rr < — — - . 

" - 2 2jt' 

Assuming that this highest axial wavenumber is well above the acoustical 
wavenumber for the highest frequency (usually a good assumption in shells), 
this gives an upper bound for fc' = « tt/Rz. Taking the common 

logarithm of both sides of 8.6 and substituting for k^, we obtain 

Rz = 207rWio(e)(^) « 27.3 (^) . (8.8) 

This expression gives us the guidance we need to choose d based on our 
knowledge of D and the shortest likely wavelength on the structure, Xmin- 
As an example, assuming D = 55 dB (about right for a 12-bit system with 
good positioning), equation 8.8 yields d « Ani,„. This could be very small 
for a thin shell at high frequencies—possibly just a few millimeters! This 
is where the process of balancing the hologram parameters for experimental 
design begins. 

How does this impact Genah ? For the most part, the quality of the 
data (its dynamic range) will influence both the maximum distance over 
which holograms may be back propagated and the magnitude of the cutoff 
wavenumber for the wavenumber filters. If the reconstruction distance is 
fixed (say, back to the source surface), then the impact is felt in the filtering 
and thus in the available resolution. If d is larger than the data quality can 
support, the Genah user may find that the filters must be so restrictive that 
the reconstructable resolution and possibly the frequency range become more 
limited than desired. This is why it is so valuable to preview the data before 
reconstruction and to perform forward-transform-only analysis (Option 1) 
to examine the wavenumber content and its noise level. 
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Chapter 9 

Analysis of FORTRAN Code 
in GENAH 


In this chapter, we examine the FORTRAN code in Genah and discuss ex¬ 
actly how it implements the holographic algorithm discussed in the previous 
chapters. For a list of subroutines and their arguments, see Appendix D. 


9.1 Variable Declarations 

Genah uses the implicit none qualifier to force explicit declaration of all 
parameters. All of the subroutines provided by NRL follow this convention. 

The first set of declarations is a group of parameters. The first four control 
array sizes: nrkcbins(250) is the number of wavenumber filters that can be 
read in; maxn(34) is the maximum order of Bessel function needed (it should 
be Timax-^^, where Umax is the highest positive circumferential order allowed); 
ldz(256) is the leading dimension of the hologram data arrays (the maximum 
number of axial points allowed); and sdz(128) is the second dimension of 
the hologram data arrays (the maximum number of circumferential points 
allowed). The other three parameters are: twopi, simply the value of 27r; rtd, 
a radians-to-degrees conversion factor; and rhoc, the specific impedance of 
water in metric units. 

Following the parameter declarations are the declarations for all of the 
internal variables used by Genah grouped in order of size; Integer, Real, 
Complex, Real*8 and Complex*16. Within these groups, variables are de- 
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dared roughly in alphabetical order, but also grouped according to usage 
in the program. The data arrays are all complex. Double precision is used 
only in computing the Hankel functions—for large orders and arguments, the 
Bessel functions grow quite large, but their ratio (needed for the propagator, 
see Section 7.4) remains of small order. 

The single INCLUDE file for Genah is the set of declarations and common 
blocks for the data file header. They are numerous and confusing. For 
clarity’s sake, they were collected into their own file. For more information 
on the Header and its contents, see Appendix B. 

The variables used in the Radix-2 Fast Fourier Transform (FFT) in 
Genah are gathered separately from the general declarations for ease of ma¬ 
nipulation. Note that while the arrays are of length Idz, the indexing vector 
i2 is of length log 2 {ldz). Since the compiler cannot perform the logarithmic 
operation, this value must be altered by hand if Idz is changed. 

Finally, the two functions cang and windo are declared external to satisfy 
the explicit declaration criterion. 


9.2 Input Data File and Parameters 

The next section of Genah concerns opening the input data file (pressure 
holograms) and setting parameters, many of which are flags which control 
decisions in the processing operations. The details of the inputs expected 
from the user axe not covered here; they axe discussed at length in Appendix 
C. However, it is appropriate to discuss some of the variables and their 
settings here. 

There are a set of counters which describe the extent of the hologram 
mesh to the various DO-loops throughout the code. In particular, neixpts, 
the number of points in the axial direction (also the leading dimension of 
the data arrays), has three versions. This is because, as a result of the zero 
padding option, the input, working and output arrays may all have different 
lengths. The second dimension of the array, ncirpts does not vary in this 
way. 

The choice of axial spatial windowing determines how the array tuk is 
loaded. “No window” actually creates the so-called boxcar function, with 
amplitude 1.0 over the aperture. The Tukey, Hanning and Hamming windows 
are discussed in more detail in Section 5.2.1. 
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Depending on the choice for zero padding, the number of axial points 
to be used in processing is set. Given that and the physical lattice spacing 
alat, the wavenumber increment delk is set according to the relation AA:^ = 
27r/(iV^ai). 

The frequency range and increment available in the hologram data file are 
announced. Based on these, the user can request which holograms should be 
processed. The resulting inputs become the frequency loop parameters. 

The next section of code deals with setting the wavenumber domain fil¬ 
ter parameters, discussed in Appendix E. If the option for a frequency- 
dependent filter file is chosen but the number of filters in the file exceeds 
the array limits set by parameter nrkcbins, a fatal error is flagged and the 
program is terminated. Also, circumferential order “low-pass” filtering can 
be set here. 

The next section loads the all-important flag if qout, the process control 
flag. The user will note that it is used throughout the remainder of the 
Gen AH code to decide where the processing path should go next. While 
the current implementation of this decision tree is archaic, it has proven 
to be robust. It may be altered in a future release. Based on the values 
of if qout and the zero padding flag izadd, the output leading dimension 
counter naxpts.out is set. 

Next, the physical parameters which describe the hologram and the re¬ 
construction surface are listed, and the user can alter them if necessary, for 
example if the reconstruction surface is to be at some radius other than the 
source surface. Again, more details axe available in Appendix C. 

When the user gives the name for the output data file(s), the code formats 
the names with the appropriate extension. It then opens the files for writing. 
Note that the open status is set to ‘unknown’: if the file exists, it will be 
overwritten; if it does not exist, it will be created. The header block is then 
loaded with the input hologram file header data after modifying it appropri¬ 
ately. To avoid conflict with common block quantities (such as nsteps and 
f start), temporary variables are used to hold the input file values until the 
output file headers have been written. 
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9.3 The Main Loop 

This is the heart of the Genah program; the loop over the individual fre¬ 
quency bins for the processing of each hologram. The loop (and thus the 
current frequency) is controlled by the parameters ihstaxt (the index of the 
first hologram of those available to be processed), ihno (the total number of 
holograms to be processed) and ihstep (the step interval, or the increment 
to the hologram index variable ihol). Note that the top hologram is not 
ihno itself, but (ihno - 1) * ihstep + ihstart. 

9.3.1 Filter Pzirameters 

If the wavenumber filter is constant (independent of frequency), it is set. 
Otherwhise, the available list of filter parameters is searched for those which 
apply to the current hologram. This is done by finding the two bin pointers in 
the filter file which bracket the current hologram index. After the hologram 
index reaches the top bin pointer, it is no longer reset; the filter parameters 
remain set at the values specified for the top bin pointer. 

9.3.2 Current Frequency and Wavelength 

The current hologram index is used to compute the center frequency of the 
hologram frequency bin, based on the frequency of the first hologram in 
the input file and the bin width of that file as freq = delf_in * (ihol 
- 1) + f St art-in. The wavelength of an acoustic wave in water at that 
frequency is then computed using a sound speed of 1480.82 m/s (clean water 
at atmospheric pressure at 20 degrees Celsius). If the data from the input 
file is taken under other conditions, the value of the sound speed should be 
altered. It is conceivable that an open slot in the header (see Appendix 
B could be used to store the temperature and depth of water at the time 
the data was taken. Since these are closely controlled in NRL’s laboratory 
environment, they have been held constant. 

9.3.3 Reading in the Hologram 

At this point, the data for the hologram pointed to by index ihol is loaded. 
The outermost IF block allows for the number of axial points to be any 
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power of 2 up to and including 256. If other values are needed, this section 
must be modified. The first step within each block is to call the subroutine 
filread_ieee. 

For any number of axial points M less than 256 (but still a power of 2), a 
temporary array zdummy is filled with the input data. If the zero padding flag 
izadd is set, the data array zdat is padded with M/2 zeros at the beginning 
and end of each scan line, and the data, multiplied by the windowing vector 
tuk, is loaded into the center M values. Otherwise, the data is windowed 
and loaded directly into the first M locations (for each column) of zdat. 

The case of 256 axial points is a bit different. Since the data arrays are 
limited to 256 values in the leading dimension, scan lines of this length cannot 
be zero padded. Also, the temporary array zdummy is not used. Data are 
loaded directly into zdat and then windowed with vector tuk. 

9.3.4 Conjugation, Maximum Search and Wavenum¬ 
ber 

With the data loaded and windowed, the hologram is searched for the data 
point with maximum modulus, the index of its location in the mesh, and 
the complex phase of that value. This information is announced to the user’s 
screen and copied into the ASCII output information file for future reference. 

If the data need to be conjugated (see Section C), it is done in this loop. 
Use of the IF test on the conjugation flag iconj inside the loop is inefficient, 
but since this loop is executed only once per hologram, it does not matter 
very much. 

Before the hologram processing begins, the current fluid wavenumber is 
computed, as well as some constants based on it that come into play in the 
propagation section. Also, the source and hologram radii are reset here. 


9.3.5 Forward FFT, Axial Direction 

The data must now be forward transformed to the wavenumber domain. To 
reduce the time spent in the transform process to a minimum, all of the 
Fourier weights (cosine and sine factors) and the butterfly bit reversal index 
are precomputed in routine wtrev. The “direction switch” isw is used to 
mark the sense of the complex coefficients and thus whether it is a forward 
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or backward transform. In this case, setting isw to 0 indicates that it is a 
forward transform. 

NOTE: The sense of the forward transform for spatial data is held to be 
the opposite of that for temporal data. This is a critical issue. An outward 
travelling wave is described by the function exp\—i{u)t — kr)]] the signs of the 
transform kernels are opposite for temporal and spatial data. This 
must be taken into account in preprocessing data for Genah ! 

The forward transforms (and all subsequent transforms) are performed 
on all the columns and then all the rows, taking advantage of the linearity of 
the Fourier transform operator in computing the two dimensional transform. 
The axial-coordinate data corresponds to the columns of the data array, and 
they are transformed first. Since the FFT routine employed here is in-place 
(destroying the input vector and replacing it with the result), the columns 
(rows) are loaded one at a time into a transform vector zf. The vector is 
transformed and loaded back into a result array. 

In loading the result array, the vector is rotated by half the aperture, 
and all the even indices are multiplied by -1. The first process rotates the 
wavenumber domain so that the zero wavenumber is in the center of the 
array rather than in the corner. The second process, which is equivalent 
to multiplication by the Nyquist wavenumber, is then needed to properly 
phase all the contributions for plotting purposes. Otherwise, the signs of 
the components toggle positive and negative. This is self-correcting when a 
forward transform is followed by an inverse transform, but when processing 
(such as Options 1, 2 and 4) stops after the forward transforms, the toggling 
must be removed to properly visualize the transforms. If the user were to 
always do full reconstructions (back to the spatial domains), these loops could 
be removed. 

9.3.6 Forward FFT, Circumferential Direction 

If Processing Option 4 is in effect, no forward transforms are performed in 
the circumferential direction. The code jumps to markers 79, 210 and 215, 
at which point the Pik^, 4>\ th) data are written to the output file. Likewise, 
if only a single scan line is available (ncirpts = 1), no circumferential trans¬ 
forms are performed. The code jumps to marker 250 to begin propagation 
(assumed to be for order n = 0 only). 

Otherwise, the forward transforms are performed on data in the circum- 
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ferential direction. Each row of the data array is loaded into the transform 
vector zf, transformed in place, and returned to the data array after hav¬ 
ing been rotated and toggled. (See Section 9.3.5 for an explanation of this 
process.) 

When the forward transforms are complete, the weighting and bit reversal 
vectors are reset for inverse transforms in the axial direction, using isw = 1. 
This section of code should really precede the backward transforms locally; 
it is located at this point for historical reasons. 

9.3.7 Computation of the Propagator 

If Processing Options 1 or 4 are in effect, no propagations are applied to the 
data. The code jumps to markers 210 and 215, at which point the Pn{kz\rn) 
or the P{kz,(f); rjj) data are written to the output file. 

Otherwise, this is the section in which the propagators for the current 
hologram frequency are generated and their dot product with the wavenum¬ 
ber domain data is formed. Propagators are computed as vectors whose 
indices are the circumferential orders n. The outer loop is over axial wave- 
number index. The axial wavenumber is computed using this index and the 
interval wavenumber delk. The square of the radial wavenumber, called tst 
here, is computed. It serves as the switch for propagating or evanescent 
waves (see Section 7.4). 

The first block executes if the radial wavenumber is real (propagating 
waves). In that case, the (VAXMATH) Bessel functions dbesy and dbesj are 
called to load vectors by, byO, bj and bjO at the hologram and source radii, 
respectively. The functions are computed up to order maxn, which is a pa¬ 
rameter set to Inwall + 2 (see Section 9.1). The Hankel function vectors are 
then formed according to Hn{krr) = JniKr) + iYn{kj.T) and stored in bh and 
bhO. 

An inner loop over circumferential orders is executed next. A pointer 
nn into the Hankel function vectors is created to translate between the loop 
counter and the order comprising that row of the data array. (Recall that the 
propagators are symmetrical in wavenumber about the zero wavenumber.) 
The circumferential wavenumber is computed for the sake of the wavenumber 
filter. An offset is added to the Hankel pointer since FORTRAN arrays are 
numbered from 1, not 0. The pressure propagator grnp is computed as the 
product of the ratio of Hankel functions (see Section 7.4) and the wavenumber 
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filter windo for the current axial wavenumber and circumferential order (see 
Section E). The propagated pressure value for the current circumferential 
order and axial wavenumber is then computed. 

The velocity propagator grnv is computed in a similar manner, except 
that the formula is slightly different for n = 0 than for n ^ 0. This is 
due to the fact that the numerator of the propagator is the first derivative 
of the Hankel function. This is most easily computed using the recursion 
relation H'^{z) = (l/2)(£r„_i(z) - Hn+i{z)) for n > 0 and Hq{z) = -Hi{z). 
The reader will also note the presence of the constant ikr/2pck^ satisfying 
Newton’s law for conversion from pressure to velocity (see Section 7.2.2). 
Again, the function windo is called to compute the value of wavenumber 
filter for the current wavenumbers. The reconstructed velocity value is then 
computed. 

Because the propagator is symmetric, the outer loop counts over only half 
the axial wavenumbers. In this inner loop, the symmetric axial wavenumber 
component is also computed using the current values of the propagators. 

The second block of this section is used if the radial wavenumber is imag¬ 
inary (rk < 0) and the components are evanescent. In this case, calls are 
made to dbesk to compute the Bessel functions of the Third Kind. Other¬ 
wise, the block proceeds identically to that for the propagating waves as just 
described. 

9.3.8 Backward FFT, Axial Direction 

If Processing Option 2 is in effect, no inverse transforms are performed. The 
code jumps to marker 215, at which point the K(A:^, n; rs) are written to the 
output file. 

Otherwise, the inverse transforms are begun on data, first in the axial 
direction. At this point, both pressure (zp) and velocity (zv) data must 
be transformed, so there are two sets of transforms within the loop on cir¬ 
cumferential order. Within that loop, if Processing Option 5 is in effect, no 
inverse transforms are necessary on pressure, since only velocity is output. 
The code jumps to marker 230, at which point the inverse transforms on 
velocity proceed. The transforms themselves are identical in structure to the 
forward transforms, discussed above, with the exception that the direction 
switch isw is set to 1 for inverse transforms. 
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9.3.9 Backwzird FFT, Circumferential Direction 

Following the inverse transforms in the axial direction, the code jumps to 
marker 215 for Processing Option 3, at which point the V{z,n\rs) data are 
written to the output file. Likewise, if only one scan line is being processed 
(synthetic axisymmetry), the code jumps to marker 250. At this point inten¬ 
sity and power are computed, and pressure, velocity and intensity are written 
to their respective output files. 

The weighting and bit-reversal vectors are next loaded for the inverse 
transform on circumferential orders. As in the axial inverse transform, both 
pressure and velocity must be treated, so there are two sets of transforms. If 
Processing Option 5 is in effect, no inverse transforms are necessary in the 
circumferential direction. The code jumps to marker 240, at which point the 
individual orders are written to the output file. This is discussed in more 
detail in Section 9.3.10 below. 

It is not immediately obvious, but only Processing Option 0 (full recon¬ 
struction) has survived to the point of inverse transforms on circumferential 
order. The transforms proceed exactly as they had in the forward transforms 
on 4>, except that the direction switch isw is set to 1. 

When the transforms are complete, the code jumps to marker 250, at 
which point the intensity and power are computed and pressure, velocity 
and intensity are written to their respective output files. This is described 
in detail in Section 9.3.11 below. 

9.3.10 Limited Order Velocity Reconstruction 

Processing Option 5 differs markedly from the other options. Instead of writ¬ 
ing a single reconstruction for each frequency, this option writes nmaix — main -I-1 
velocity reconstructions for each frequency. However, the reconstruction is 
not a function of circumferential order itself, that is, V{z,n;rs) as in Op¬ 
tion 3. Instead, for each order nn between main and nmax, both the positive 
and negative orders are loaded by themselves into an otherwise empty vector 
and inverse transformed. This yields the reconstructed velocity Vn{z, <t>'\ rs), 
where the (j)' is used to indicate that only orders ±n are present in the velocity 
field. 

Again, it is emphasized that for each frequency bin. Processing Option 5 
writes nmax — main -I-1 entries into the output file. Thus, output files will 
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be much larger in size than would otherwise be the case. More importantly, 
when plotting or post-processing, the user should remember that only the 
first reconstruction will have the correct frequency associated with it. The 
number of orders for each bin would have to be used to compute its frequency, 
which is not a standard procedure in most post-processing programs. 

When the inverse transforms are completed for each order of interest, the 
reconstructions are written to the output file. The code which does this is 
discussed in more detail in Section 9.3.12 below. 


9.3.11 Intensity and Radiated Power 

The acoustic intensity is computed exactly as prescribed in equation 7.36 
of Section 7.5.2. At each lattice location (i, j), the value of the pressure 
zp(i,j) is multiplied by the conjugate of the value of the radial velocity 
zv(i, j) and reduced by half. The complex result is stored in array zdat. 

Prior to the double loop over the aperture, the accumulator power and a 
maximum intensity store pmax are zeroed. Within the loops, the real part of 
the intensity at each lattice point is extracted from zdat(i, j) and stored 
in ptenn. This term is compared against pmax; if it exceeds pmaix, the value 
of pmax is reset to pterm, the current (now maximum) intensity. The value 
of pterm is then added to the power accumulator. Note that this is a crude 
rectangular quadrature approximation to the integral of equation 7.37. How¬ 
ever, as discussed in Chapter 8, the mesh should be fine enough to adequately 
represent the wave structure of the field. This implies sampling sufficient to 
allow the rectangular quadrature to approximate the power integral. After 
the loops are complete, the power accumulator is multiplied by Az (alat) 
and A0 = 2T:r/n^ (radius and ncirpts). These represent the differentials 
dz and d(j) in the integral of equation 7.37. Note that both alat and radius 
are multiplied by conversion-factor: this guarantees that the power is in 
SI units regardless of the units of the distance measures. 

The maximum intensity and total power in the hologram are reported 
to the ASCII output file if Processing Option 0 (Normal Reconstruction) 
is in effect. Next, the intensities throughout the hologram are searched for 
any whose real part is within a fraction of the maximum. This fraction is 
preset in the code as parameter cutof. In the current version, cutof is 
set to 0.7 or 70% of the maximum real part of the intensity. Every point 
in the hologram whose intensity exceeds cutof has its pressure (magnitude 
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and phase), velocity (magnitude and phase), intensity (magnitude) and the 
difference between the pressure and velocity phases written to the output file 
along with the indices of its location in the mesh. The phase difference is 
useful in determining how reactive the field is; as the difference approaches 
±90°), the intensity is becoming purely reactive and no power is radiating. 
These values give an indication of how energy is distributed throughout the 
wave field. If only a few points are listed, then the maximum point (usually 
the drive point) stands alone. If many points are listed, then the source is 
probably near some sort of resonance where most of the surface is strongly 
contributing to the field. In essence, the quantities written to the ASCII 
output file give a qualitative feel for the character of the wave field. 

9.3.12 Writing the Reconstructions: Spatial Domain 

With all computations completed for Processing Option 0, the reconstruc¬ 
tions can now be written to disk. Holograms are written using subroutine 
f ilwrite.ieee. They are written in fixed record length binary format (see 
Appendix A for details). The record lengths for all holography files are fixed 
at 128 longwords (512 bytes), so different numbers of records must be used 
for different numbers of axial points. A 128-longword record will hold 64 
complex values. For cases where the number of axial points per scan line is 
less than 64 (but still a power of 2), more than one scan line can be written 
to a record, so the number of records needed to write an entire hologram is 
computed as 

num_recs = ncirpts/(64/naxpts_out) , (9.1) 

where integer arithmetic is enforced. For the cases where the number of axial 
points per scan line is 64 or 128, the number of records is simply the number 
of scan lines, ncirpts. In that case, the subroutine uses the header value 
nsteps(l), stored in the COMMON block, to determine how many records are 
needed. 

If the input data have been zero-padded, there is no need to write out the 
padded region of the reconstructions. In that case, only the interior values 
are written out. This is accomplished by loading the values into the array 
zdummy and writing that array rather than the data array. In either case, 
three arrays are written: pressure to logical unit 8, velocity to logical unit 2 
and intensity to logical unit 3. Once the data for the current reconstructions 
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are written, the code jumps to marker 220, after which the main loop over 
frequency bin is terminated. If the current frequency was the last to be 
processed, all files are closed and Genah terminates. 

9.3.13 Writing the Reconstructions: Wavenumber Do¬ 
main 

Marker 215 is the jump point for output for Processing Options 1, 2, 3 or 4 
(wavenumber domain output). As with spatial domain data, holograms are 
written using subroutine f ilwrite.ieee. The same 128-longword records 
are used for this data. Output for wavenumber domain data is the same 
as that for spatial domain data with one exception. The case of 256 axial 
points is allowed, since zero-padded data has interpolated values throughout 
the entire wavenumber domain. In other words, no distinction is made as 
to whether the data have been padded or not in the wavenumber domain. 
Processing Options 1, 2 and 4 fall into this category. In each case, the data 
are written to logical unit 2. 

The case for Processing Option 3 {V{z,n)) is similar to that for spatial 
domain data in that if the data are zero-padded, only the center values are 
written. As such, it has its own processing block. In this case, the velocity 
data are written to logical unit 2. 

Once the data for the current reconstruction are written, the code passes 
marker 220, after which the main loop over frequency bin is terminated. If 
the current frequency was the last to be processed, all files are closed and 
Genah terminates. 
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Appendix A 

The Holography File Format 


The input and output files for Genah are described generally in Section 
6. This appendix formalizes that discussion in sufficient detail to allow 
the user to create or read files in that format. If the user wishes to use 
some other scheme, the Genah code can be altered. For input, the subrou¬ 
tine filread-ieee would need to be replaced. For output, the subroutine 
f ilread-write would need to be replaced. 

All Genah data are stored in an unformatted, binary fashion. They are 
not transportable across different types of computers (such as Sun, VAX, PC 
and so on) or Operating systems (such as Unix, VMS, DOS, and so on). The 
data are written and read using the FORTRAN direct access method with a 
record size of 128 longwords (512 bytes). Note that a program cannot read 
the data using records whose length is an integer multiple of 128 longwords 
due to the presence of the header, which offsets the record pointer. The data 
can be read and written using the file manipulation routines contained in the 
stdio library of most C compilers. 

Holography data can be described in either a deconstructive or a con¬ 
structive way. From the top down, a file is made up of a header and all the 
holograms. A hologram is made up of scan lines. A scan line is made up of 
step locations or “points”. A point contains 1 complex datum representing 
the calibrated value of pressure (velocity, etc.) per unit input. From the 
bottom up, a scan line is made up of one or more records. A hologram is 
made up of the records for all the scan lines at one frequency. A file is made 
up of the records for the holograms for a set of frequencies and the header. 

A Genah record is 512 bytes long. This is declared in the FORTRAN-77 
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OPEN statement as recordsize=128 longwords at 4 bytes per longword. This 
record can represent 64 complex numbers at 2 longwords per number. Both 
input and output files are stored in direct access mode with this fixed record 
length. To access a record, a pointer (called an associatevariable in 
FORTRAN-77 ) must be declared in the OPEN statement. The OPEN statement 
should also include the keywords form='unformatted’, access= ‘direct ’. 
An example is shown in the code fragment at the end of this appendix. 

For scan lines of M < 64 steps (but still a power of 2), N scan lines are 
packed into each record according to the rule N = QifM. For scan lines of 
M = 64 steps, each scan line constitutes a single data record. For scan lines 
of M > 64 steps (but still a power of 2), AT records are used for each scan line 
where N = M/64. So, for instance, two 32-step scans would fit into a single 
record, a 128-step scan would require two records and a 256-step scan would 
require four records. A 154-step scan would require three records, although 
the third record would be padded out with 38 empty entries. The number 
of records in each hologram is equal to the number of scan lines times the 
number of records per scan line. A typical cylindrical hologram of 128 points 
per scan line and 64 scan lines would require 128 records of 512 bytes each, 
for a total of 64 kilobytes of storage. 

Each file written in NAH format has a 512-byte header. Record number 
1 of each file (whether input or output) is dedicated to this header record. 
The contents of the header are described in detail in Appendix B. 

The size of a file in bytes is exactly computable. It is equal to the number 
of holograms (frequency bins) times the number of scan lines per hologram 
times the number of records per scan line times 512 bytes per record plus 
512 bytes for the header. 

The following code fragment demonstrates the key elements in writing 
data to a file in Genah format. Reading data is exactly the same, where 
write statements are replaced by read statements (see the Genah subrou¬ 
tine filread-ieee.f. 

The code is self-explanatory with the exception of variables il and i2. 
It is up to the user (code) to decide which values along the first coordinate 
of the hologram should be written at each loop. For example, suppose a 
cylindrical hologram has 128 steps per scan line (two records for each line) 
and all values should be written. The loop on k will execute twice. On the 
first pass, first_x(l) should generate the value “1” and last_x(l) should 
generate the value “64”. On the second pass, first_x(2) should generate 
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the value “65” and last_x(2) should generate the value “128”. Suppose, 
however, that the internal representation of the hologram has been padded 
out to 256 steps per scan line, but only the center 128 values occur over 
the original hologram aperture. Only two records are required for the 128 
points of interest, but first_x(k) should generate the values “65” and “129”, 
respectively, and last_x(k) should generate the values “128” and “192”, 
respectively. Then the control values il, i2 will properly point into the 
hologram array in the write statement’s implied DO loop. Note that the 
current version of Genah applies a brute force method of accomplishing this 
based on the processing path. A more flexible method will be used in an 
upcoming release. Note, too, that this discussion only applies to writing a 
hologram; when reading, the hologram array is just loaded beginning at the 
first value along scan lines. 


c 

c ” Open the file as direct access, unformatted, 512-byte records, 
c 

open (unit=iunit, file=filename, status=’unknown’, 

& form=’unformatted’, access=’direct’, recordsi 2 e= 128 , 

& associatevariable=nrec) 

. (Load parameter values into 512-byte header array) 


write (iunit, rec=l) ihead !Write the header into record #1 

do hol_count = 1, nholograms !Loop over frequency bins 
. (Process holograms) 


c 

c 

c 


do j = 1, nscan.lines 

do k = 1, recs_per_scan 

11 = first_x(k) 

12 = last_x(k) 


ILoop over scan lines 
ILoop over multiple records/scaii line 
!User must decide which points get 
! written from a scan line. 


Set the record number pointer. 


nrec = ((hol.count - 1) * recs_per_scan * npts_y) 
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c 

c 

c 


& + recs_per_scan * j 

& + (k - 1) 

— Write the record using Implied DO along the 1st coordinate. 


write (iunit, rec=nrec) (zhol(i, j), i = il, i2) 


end do 
end do 
end do 


!End loop over records/scaji line 
!End loop over sccin lines 
!End loop over frequency bins 


92 



Appendix B 

The Holography File Header 


B.l Introduction 

This section describes in detail the structure and contents of the header block 
prepended to any Genah data file. The block is one record long, for a total 
of 512 bytes. (See Appendix A for more information about records.) 

The header was originally developed for use on a PDP-11 series computer 
whose native word length was 16 bits, not 32 bits as in most workstations 
in use today. As a result, integers were stored with only 2 bytes. It is not 
difficult to handle the conversion to longwords (4 byte integers) in FORTRAN 
. There is an unfortunate side effect to the use of the shorter integers: most 
of the real numbers (stored in 4 bytes as they are on 32-bit machines) do 
not begin on a longword (4-byte) boundary. While VAX systems will au¬ 
tomatically compensate for this (thus providing backward compatibility to 
PDP-11’s) through the FORTRAN EQUIVALENCE statement, most systems will 
not. While many vendors, such as SGI and SUN, provide a library of calls to 
handle misaligned data, the traps and handlers they use are very expensive 
at runtime. As a result, the read and write routines in Genah simply use 
a BYTE array of length 512 for the header block. The parameters within the 
header are then individually equated with specific locations in the array, and 
the pass-by-reference standard of FORTRAN-77 handles the values automat¬ 
ically. While this makes for lengthy coding in those subroutines, it works 
quickly and efficiently at runtime. 

All of the parameters stored in the Genah header are listed below in 
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Table B.l. This is followed by an in-depth discussion of each parameter and 
its use. 


B.2 Description of Header Parameters 

The following is a detailed description of each of the parameters currently de¬ 
fined in the Genah file header. Each parameter is itemized by its FORTRAN 
name as it appears in the Genah code. The parameter’s function is de¬ 
scribed, followed by the range of allowable values, the number of bytes (in 
decimal) the variable is offset in the 512-byte header, and the FORTRAN vari¬ 
able type used to represent the parameter. 

Parameters were listed strictly in offset order in the previous Table B.l 
for reference purposes. Note that they are grouped here by functionality, 
not necessarily by order of their offset in the header. While the rough order 
of their appearance in the header is followed, the user will see deviations 
from that order when a parameter’s function is more closely related to non¬ 
neighboring parameters. 

B.2.1 Experiment, Target and File Description 

These entries in the file header are used to provide an accurate description 
of the experimental run, the target under test, and the state of the data 
contained within the file. The descriptive strings can be whatever the user 
chooses; formation of an in-house standard for naming and numbering exper¬ 
iments and targets is recommended. The date and time entries refer to the 
file creation, not the execution of the experiment. They take the standard 
numerical forms MM-DD-YY and HH:MM:SS, although the date is actually nu¬ 
meric (integers), while the time is ASCII-encoded characters. This conforms 
to the FORTRAN-77 system TIME and DATE utilities. 

The radius of the target shell is included here since it forms a critical part 
of the description of the target. 

The flag NFLAG is used at NRL to describe the status of the data. In 
combination with the file name extension, it indicates what processing has 
been completed on the data up to this file. It is purely an in-house standard; 
the user may use any code system desired. The codes used by NRL are 
described in Table B.2.1. One alternative suggestion would be to use the 
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Byte Offset 


000 

050 

052 

054 

056 

058 

060 

062 

100 

102 

104 

106 

112 

114 

116 

118 

120 

122 

124 

174 

176 

178 

180 

182 

188 



Data Type I Name Variable Description 


TITLE Name, description of run 

MONTH 2-digit Month Code (01-12) 

DAY 2-digit Day Code (01-31) 

integer*2 YEAR 2-digit Year Code (83-92) 

character*2 HOUR 2-char. Hour Code (00-24) 

character*2 MINUTE 2-char. MinuteCode (00-59) 

character*2 SECOND 2-char. SecondCode (00-59) 

character*38 TARGET Name and description of target 
integer*2 NFQC *Center Frequency mantissa 

integer*2 NFQCXP ^Center Frequency exponent 

integer*2 NUMREC # of Data Records in file 

character*6 DATATYPE 6-character descriptor of data type 

integer*2 NDCVAL *# Bits DC Offset from calibration 

NCKRT *Digitizer Rate mantissa 

NEXP ^Digitizer Rate exponent 

NAVGS *# of records averaged per sample 

integer*2 NWVSCL *Peak Output Drive Voltage xlO 

integer*2 IRATE *Drive Signal Output rate (us) 

character*50 SOURCE *Name of Drive Waveform file 

integer*2 NXO Sensor start position, x Coord. xlOO 

integer*2 NYO *Sensor start position, y Coord. xlOO 

integer*2 NZO *Sensor start position, z Coord. xlOO 

integer*2 NTO *Sensor start position, 0 Coord. xlOO 

integer*2 NSCANS(3) # Scan surfaces in x, y, z or z^ 0^ r 

integer*2 NSTEPS(3) # Steps per scan in x, y, 2 or z, 6, r 



Table B.l: NAH Data File Header Block Layout. 

{Note 1: The header is a total of 512 bytes in length.) 
{Note 2: A* indicates that the field is not always filled.) 
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Byte Offset 

Data Type 

Name 

Variable Description 

194 

integer*2 

lAXISO) 

*Scan Axis Flags (Step, Scan, Surf.) 

200 

integer*2 

IDIR(3) 

*Start Direction Flags 

206 

real*4 

DISKS) 

Step, Scan and Surface step sizes 

218 

character* 18 

STRING 

*3 5-char, settling times between steps 

236 

character* 20 

MOVEFILE 

*Name of Motion Control parameters file 

256 

integer*2 

KBITS 

*# of Bits in each Digitizer Channel 

258 

integer*2 

NPOINTS 

*# of Data points per Record 

260 

integer*2 

NANGLE 

* 7 ^ of Angles in farfield recons. 

264 

268 

(6 Empty Bytes) 
integer*2 

NSCNTP 

*Scan Type Code 

270 

real*4 

GAINS(2) 

*Total Channel Gain 

278 

real*4 

SNSTVT(2) 

*Transducer Sensitivity 

286 

real*4 

XSTART3D 

*Origin value of 1st Data Coordinate 

290 

real*4 

DX3D 

*Interval value of 1st Data Coordinate 

294 

real*4 

YSTART3D 

*Origin value of 2nd Data Coordinate 

298 

real*4 

DY3D 

*Interval value of 2nd Data Coordinate 

302 

real*4 

ZSTART3D 

*Origin value of 3rd Data Coordinate 

306 

real*4 

DZ3D 

*Interval value of 3rd Data Coordinate 

310 

real*4 

FSTART 

Frequency of first Bin 

314 

real *4 

FDIF 

Frequency interval between bins 

318 

real*4 

RADIUS 

Radius of target surface 

322 

real*4 

DELTHETA 

*Polar step angle in farfield recons. 

326 

real*4 

DATAMAX 

*Largest data value (mag.) in file 

330 

integer*2 

IBINDATAMAX 

*3rd Coord. Index of DATAMAX 

332 

integer*2 

IXDATAMAX 

*lst Coord. Index of DATAMAX 

334 

integer*2 

lYDATAMAX 

*2nd Coord. Index of DATAMAX 

336 

508 

(172 Empty Bytes) 
integer*2 

NBINS 

Number of 3rd Coord. Intervals 

510 

integer*2 

NFLAG 

3-Digit Process Type Code 


NAH Data File Header Block Layout, continued. 
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available 16 bits as individual flags indicating processing which has taken 
place so far on the data. For example, the first bit could indicate time 
or frequency domain, the second and third bits could indicate spatial or 
wavenumber domains, and so forth. 

TITLE String containing the user-created name of an experiment run and 
possibly some description of it. 

Allowable Values: Up to 50 alphanumeric (ASCII) characters 

Byte offset: 000 

Variable Type: Character*50 

Criticality: Useful but not critical 

MONTH 2-digit code representing the Month portion of the date and time the 
file was created. 

Allowable Values: 01 through 12 

Byte offset: 050 

Variable Type: Integer*2 

Criticality: Automatically generated by Genah 

DAY 2-digit code representing the Day portion of the date and time the file 
was created. 

Allowable Values: 01 through 31 
Byte offset: 052 
Variable Type: Integer*2 

Criticality: Automatically generated by Genah 

YEAR 2-digit code representing the Year portion of the date and time the file 
was created. 

Allowable Values: 00 through 99 

Byte offset: 054 Variable Type: Integer*2 

Criticality: Automatically generated by Genah 

HOUR 2-character code representing the Hour portion of the date and time 
the file was created. 

Allowable Values: ASCII characters 00 through 24 

Byte offset: 056 

Variable Type: Character*2 

Criticality: Automatically generated by Genah 


97 




MINUTE 2-character code representing the Minute portion of the date and 
time the file was created. 

Allowable Values: ASCII characters 00 through 59 

Byte offset: 058 

Variable Type: Character*2 

Criticality: Automatically generated by Genah 

SECOND 2-character code representing the Second portion of the date and 
time the file was created. 

Allowable Values: ASCII characters 00 through 59 

Byte offset: 060 

Variable Type: Character*2 

Criticality: Automatically generated by Genah 

TARGET String containing the user-created name of the target under study 
in current experiment and possibly a description of it. 

Allowable Values: Up to 38 alphanumeric (ASCII) characters 

Byte offset: 062 

Variable Type: Character*38 

Criticality: Useful but not critical 

RADIUS Radius of cylindrical target’s outer surface measured from the axis 
of symmetry in the current engineering units. 

Allowable Values: Single-precision fioating point value (positive) 

Byte offset: 318 
Variable Type: Real*4 
Criticality: Critical 

NFLAG 3-digit code representing the type of processing used to create data 
in this file. Allowable Values: (See Table B.2.1) 

Byte offset: 510 
Variable Type: Integer*2 

Criticality: Automatically generated by Genah 


B.2.2 Number and Type of Data Records 

These parameters describe the size of the file to the calling program. With 
knowledge of the number of records, the data type and the number of points 
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Flag 

Value 

Processing type 

Data type 

Processing 

Generation 

444 

Structural Intensity 

Complexes 

4 

445 

Structural Power 

Complexes 

4 

554 

Acoustic Intensity I{z,((>-,r') 

Complexes 

3 

555 

Velocity V{z,4)]r'), V{z,4)'‘,r') 

Complex* 8 

3 

556 

Pressure P{z, 0; r'), P{k 2 , n; vq) 

Complex*8 

3 

557 

Velocity V(k^, n; r') 

Complex*8 

3 

558 

Velocity V{z, n; r') 

Complex* 8 

3 

559 

Pressure p{k 2 ,4>\ tq) 

Complex*8 

3 

666 

Pressure P(z,(j)-,rQ) 

Complex*8 

2 

776 

Force Reference . REF 

Complex* 8 

1 

777 

Calibrated Pressure . CPX 

Complex* 8 

1 

887 

Raw A/D, Pressure .PRS 

Integer*2 

0 

888 

Raw A/D, 2 Channels .SCN 

Integer*2 

0 

99x 

Scattered Holograms 

Complex*8 

3 


Table B.2: Values for processing type code NFLAG. 


per record, the routine knows the size of the file. In fact, this is not often 
necessary for follow-on processing, so these parameters can be considered 
optional. The parameter NPOINTS is usually used to store the number of 
time samples per record in the original acquisition of hologram data. This 
can be used along with the sample rate to determine the frequency bin width 
and available frequency range. 

NUMREC Total number of data records in file. 

Allowable Values: 1 through 32767 
Byte offset: 104 
Variable Type: Integer*2 
Criticality: Critical 

DATATYPE 6-character descriptor of data word type. 

Allowable Values: CMPLX for Complexes data, REAL for Real*4 data, 

INTGR4 for Integer*4 or INTGR2 for Integer*2 

Byte offset: 106 

Variable Type: Character*6 

Criticality: Useful but not critical if data type is known 
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NPOINTS Number of digitized time samples in original data. 

Allowable Values: 1 through 32767 
Byte offset: 258 
Variable Type: Integer*2 
Criticality: Critical 

B.2.3 Digital Acquisition Parameters 

With the exception of NCKRT, these parameters are used only when the data 
are not yet calibrated. They provide all of the channel information neces¬ 
sary to convert from digital samples to engineering units. If the data have 
already been calibrated (usually the case), only NCKRT is needed. It gives 
the digitization period (or rate) in microseconds; half of the inverse of this 
number represents the highest allowable frequency (seldom the highest useful 
frequency) in a data set. 

NDCVAL Number of bits offset for zero-volt input (DC) at A/D input, from 
calibration. 

Allowable Values: 0 through 32767 
Byte offset: 112 
Variable Type: Integer*2 

Criticality: Not critical if data are already calibrated 

NCKRT Mantissa of input A/D digitizer clock rate in microseconds. 
Allowable Values: 1 through 32767 
Byte offset: 114 
Variable Type: Integer*2 
Criticality: Critical 

NEXP Exponent of input A/D digitizer clock rate. 

Allowable Values: -32767 through 32767, usually -6 through -3 
Byte offset: 116 
Variable Type: Integer*2 

Criticality: No longer used; rate is assumed to be in microseconds 

NAVGS Number of time series averaged per sample. 

Allowable Values: 1 through 32767 
Byte offset: 118 
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Variable Type: Integer*2 
Criticality: Not critical 

NBITS Number of bits in each digitizer channel. 

Allowable Values: 10 through 16 
Byte offset: 256 
Variable Type: Integer*2 

Criticality: Not critical if data are already calibrated 

GAINS(2) Total channel gain from transducer to A/D, two channels. 
Allowable Values: Single-precision floating point value 
Byte offset: 270 
Variable Type: Real*4 

Criticality: Not critical if data are already calibrated 

SNSTVT(2) Transducer sensitivity in engineering units, two channels. 
Allowable Values: Single-precision floating point value 
Byte offset: 278 
Variable Type: Real*4 

Criticality: Not critical if data are already calibrated 

B.2.4 Source Waveform Information 

These parameters characterize the source signal used to drive the target or 
to insonify it. None of these values is critical to Genah since the source 
characteristics have already been convolved out of the pressure signal be¬ 
fore reconstructions are begun. They are carried in the header mostly as a 
reference for the user. 

When broadband nearfield acoustical holography was first developed, the 
drive waveform was a simple pulse, constructed to have a very broad fre¬ 
quency response. This response was the main lobe of a Kaiser-Bessel win¬ 
dowed sinusoid. As such, the lobe was well characterized by its center fre¬ 
quency. This value is contained in the two parameters NFQC and NFQCXP. 
(Note that the use of integers to store real values in mantissa-exponent form 
predates use of the VAX and addition of other, real-valued parameters to 
the header.) The center frequency of the drive waveform is computed as 
fc = NFQC X 10.Note that recent holograms have been generated 
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using chirped pulses with a nearly flat response across the band of interest. 
As a result, these parameters are seldom used anymore. 

The output drive voltage is measured after the last stage of amplification, 
that is, it is the actual voltage appearing at the driver input. The integer 
value in NWVSCL must be divided by 10 to get the value in volts. The drive 
signal output rate refers to the clock rate, in microseconds, of the D/A con¬ 
verter used by the arbitrary waveform generator. Finally, if the arbitrary 
waveform is stored in a computer file, the name of that file can appear in the 
string SOURCE. 

NFQC Mantissa of the source waveform center frequency. 

Allowable Values: 0 through 32767 
Byte offset: 100 
Variable Type: Integer*2 
Criticality: Not critical, no longer used 

NFQCXP Exponent of the source waveform center frequency. 

Allowable Values: -32767 through 32767 (usually less than 10) 

Byte offset: 102 

Variable Type: Integer*2 

Criticality: Not critical, no longer used 

NWVSCL Peak-to-peak output drive voltage xlO. 

Allowable Values: 0 through 32767 
Byte offset: 120 
Variable Type: Integer*2 
Criticality: Not critical 

IRATE Aribitrary waveform generator output clock rate (ps). 

Allowable Values: 0 through 32767 
Byte offset: 122 
Variable Type: Integer*2 
Criticality: Not critical 

SOURCE Name of file containing Drive Waveform. 

Allowable Values: Up to 50 alphanumeric (ASCII) characters 

Byte offset: 124 

Variable Type: Character*50 

Criticality: Not critical 
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B.2.5 Hologram Aperture Information 

As might be expected, a large number of parameters are used to describe 
the type, size and extent of the hologram aperture itself. Many of the pa¬ 
rameters date back to a period when many different types of apertures were 
being explored, including planar, modified planar, cylindrical and multiple 
surface scans. To accomodate all of the options considered possible at that 
time, an extensive set of parameters was included in the header so that a 
“smart program” could reconstruct the aperture if necessary. Now that the 
bulk of experiments (particularly those analyzed using Genah ) are simple 
cylindrical single-surface scans, several of the parameters have fallen into 
disuse. This is noted where it applies. 

Of all these parameters, the most important ones include NXO, which is 
used to represent the standoff distance (known as DRHX); NSTEPS, the num¬ 
ber of points in each coordinate comprising the mesh; and DIST, the point 
separation distances in each coordinate of the mesh. Without these values, 
reconstruction cannot take place. 

These parameters depend on the scans occurring in separable coordinate 
systems with equally-spaced points, as does Genah itself. For conformal 
holography scans, where the mesh bends around the target and the step sizes 
change, the parameter MOVEFILE becomes critical, because the file named by 
this string is used to hold the point coordinate locations. Again, this is not 
as important for Genah . 

NXO Sensor start position in x or radial coordinate xlOO. 

Allowable Values: -32767 through 32767 (positive when radial) 

Byte offset; 174 
Variable Type: Integer*2 

Criticality: Critical: used to store radial standoff, DRHX 

NYO Sensor start position in y or axial coordinate xlOO. 

Allowable Values: -32767 through 32767 
Byte offset: 176 
Variable Type: Integer*2 
Criticality: Not critical, seldom used 

NZO Sensor start position in ^ coordinate xlOO. 

Allowable Values: -32767 through 32767 
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Byte oIFset: 178 
Variable Type: Integer*2 
Criticality: Not critical, seldom used 

NTO Sensor start position in 0 coordinate xlOO. 

Allowable Values: -32767 through 32767 
Byte offset: 180 
Variable Type: Integer*2 
Criticality: Not critical, seldom used 

NSCANS(3) Number of scan surfaces in x, y, z or z, 6, r. 

Allowable Values: 1 through 32767, usually 1 or 2 
Byte offset: 182 
Variable Type: Integer*2 
Criticality: Not critical, seldom used 

NSTEPS(3) Number of steps per scan in x, y, z or z, 6 , r. 

Allowable Values: 1 through 32767 
Byte offset: 188 
Variable Type: Integer*2 
Criticality: First and second values critical 

IAXIS(3) Scan axis flags in (1) step direction, (2) scan line direction, (3) 
surface direction. 

Allowable Values: l = x,2 = y,3 = z,4 = 6,b = r 

Byte offset: 194 

Variable Type: Integer*2 

Criticality: Not critical if type of scan known 

IDIR(3) Start direction flags in (1) step direction, (2) scan line direction, 
(3) surface direction. 

Allowable Values: 1 in positive coordinate direction, —1 in negative 
coordinate direction 
Byte offset: 200 
Variable Type: Integer*2 

Criticality: Not critical if orientation of scan is known 

DIST(3) Distance between mesh locations in (1) step direction, (2) scan line 
direction, (3) surface direction, in engineering units. 
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Allowable Values: Single-precision floating-point number 

Byte offset: 206 

Variable Type: Real*4 

Criticality: First and second values critical 

STRING Three 5-character strings representing the time allowed for settling 
between steps in (1) step direction, (2) scan line direction, (3) surface 
direction in the format MM: SS. 

Allowable Values: 00:00 through 99:59 in ASCII numerics 

Byte offset: 218 

Variable Type: Character* 18 

Criticality: Not critical, seldom used 

MOVEFILE Name of file containing robotic manipulator’s motion control pa¬ 
rameters. 

Allowable Values: Up to 20 alphanumeric (ASCII) characters 

Byte offset: 236 

Variable Type: Character*20 

Criticality: Not critical 

NSCNTP Code indicating type of scan. 

Allowable Values: 10 for single line scan, 20 for single surface (multiple 

lines) scan, 30 for multiple surface scan 

Byte offset: 268 

Variable Type: Integer*2 

Criticality: Not critical, no longer used 

B.2.6 Farfield Reconstruction Angle Information 

These two parameters are not used by Genah . They are reserved for post¬ 
processing in which holographic data is projected to a circle in the farfield. 
They are listed here for completeness. 

NANGLE Number of angles computed. 

Allowable Values: 1 through 32767 
Byte offset: 260 
Variable Type: Integer*2 

Criticality: Farfield reconstructions only; not used in Genah 
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DELTHETA Polar angle step size. 

Allowable Values: Single-precision floating-point value 
Byte offset: 322 
Variable Type: Real*4 

Criticality: Farfield reconstructions only; not used in Genah 

B.2.7 Data Block Descriptors 

One can picture the data in a holography file as a large cube or block. Two 
of the coordinates of the cube represent the spatial dimensions spanned by 
the hologram aperture. The third coordinate represents the frequency range 
over which the holograms exist. Since all three coordinates are independent 
under Fourier Transform, each can be replaced by its transformed coordi¬ 
nate. For a program handling this “data block”, it is important to have 
information describing the coordinate origin (the “comer” of the block) and 
the sizes of the discrete steps along each coordinate. The first six parameters 
of this group fall into that category. For historical reasons, the next two 
parameters provide that information redundantly in the most common case 
of holograms at individual frequencies. FSTART is the center frequency of the 
first hologram, and FDIF is the frequency interval between bins. 

The next four parameters deal with the magnitude and location in the 
data block of the global maximum of the file. These parameters are not 
explicitly required by Genah , nor are they set by it. However, they often 
prove useful for post-processing which might include visualization, for which 
scaling to the global maximum is helpful. If it is provided in the header, then 
the (somewhat expensive) search for the maximum need only be accomplished 
once regardless of how many times the data are later processed. 

The last parameter, NBINS, is a critical one for Genah . It contains 
the number of frequency bins (holograms) present in the data file. While it 
might be considered redundant with NPOINTS, with use at NRL the former 
has come to be associated with frequency bins while the latter is associated 
with the number of time samples per record taken in the original experiment. 
The key difference is that while NPOINTS is held to the same value throughout 
processing steps, NBINS reflects the current number of frequency bins (or time 
steps in the case of back-transforming to the time domain) actually contained 
in the data file. Since the range of bins can be widely different from one file 
to the next, NBINS becomes a critical means of determining the extent of the 
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data in the current file. 


XSTART3D Origin value of first data coordinate. 

Allowable Values: Single-precision floating-point value 
Byte offset: 286 
Variable Type: Real*4 

Criticality: Not critical, not used or set by Genah 

DX3D Interval value of first data coordinate. 

Allowable Values: Single-precision floating-point value 
Byte offset: 290 
Variable Type: Real*4 

Criticality: Not critical, not used or set by Genah 

YSTART3D Origin value of second data coordinate. 

Allowable Values: Single-precision floating-point value 
Byte offset: 294 
Variable Type: Real*4 

Criticality: Not critical, not used or set by Genah 

DY3D Interval value of second data coordinate. 

Allowable Values: Single-precision floating-point value 
Byte offset: 298 
Variable Type: Real*4 

Criticality: Not critical, not used or set by Genah 

ZSTART3D Origin value of third data coordinate. 

Allowable Values: Single-precision floating-point value 
Byte offset: 302 
Variable Type: Real*4 

Criticality: Not critical, not used or set by Genah 

DZ3D Interval value of third data coordinate. 

Allowable Values: Single-precision floating-point value 
Byte offset: 306 
Variable Type: Real*4 

Criticality: Not critical, not used or set by Genah 
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FSTART Frequency of first hologram bin. 

Allowable Values: Single-precision fioating-point value 
Byte offset: 310 
Variable Type: Real*4 
Criticality: Critical 

FDIF Frequency interval between hologram bins. 

Allowable Values: Single-precision floating-point value 
Byte offset: 314 
Variable Type: Real*4 
Criticality: Critical 

DATAMAX Maximum data value (magnitude) in file. 

Allowable Values: Single-precision floating-point value 
Byte offset: 326 
Variable Type: real*4 

Criticality: Not critical, not used or set by Genah 

IBINDATAMAX Third coordinate index of DATAMAX. 

Allowable Values: 1 to NBINS 
Byte offset: 330 
Variable Type: Integer*2 

Criticality: Not critical, not used or set by Genah 

IXDATAMAX First coordinate index of DATAMAX. 

Allowable Values: 1 to NSTEPS(l) 

Byte offset: 332 
Variable Type: Integer*2 

Criticality: Not critical, not used or set by Genah 

lYDATAMAX Second coordinate index of DATAMAX. 

Allowable Values: 1 to NSTEPS(2) 

Byte offset: 334 
Variable Type: Integer*2 

Criticality: Not critical, not used or set by Genah 

NBINS Number of holograms, time steps, or third coordinate intervals. 
Allowable Values: 1 to NPOINTS 
Byte offset: 508 


108 




Variable Type: Integer*2 
Criticality: Critical 



Appendix C 

List of Inputs to GENAH 


This appendix is a compilation of Genah input parameters and allowed 
responses to them. Explanations of the parameters and their effects are not 
provided here. For more detail, a full description is available in Chapter 5. 


Parameter: iconj 
Allowable 1 
Responses: 0 


Incoming data Conjugation flag 
Conjugate all incoming data 
Do not conjugate incoming data 


Paramter: ivax 
Allowable 1 
Responses: 0 


Binary number format flag 

Convert from DEC format to IEEE format. 

Do not perform any number format conversions. 


Parameter: fname 
Allowable (Unix) 
Responses: 


(VMS) 


Input Filename 

name or 

./name or 

dir/name or 

dir/subdir/.. ./name or 

/dir/subdir/.. ./name or 

/partition/dir/subdir/.. ./name etc. 

name.ext or 

[dir]name.ext or 

[dir.subdir...]name.ext or 

disk: [dir. subdir...] name. ext or 
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machine::disk: [dir.subdir...]name.ext etc. 
(String, < 80 characters) 


Parameter: units 
Allowable 0 
Responses: 1 
2 


Type of Units for data in input file 
Inches (English) 

Meters (MKS) 

Centimeters (cgs) 


Parameter: ncirpts 

Allowable “n” 
Responses: 


Number of scan lines per hologram 
(Prompted only when header value is not 
an even number.) 

Actual number of scan lines in file 
(EVEN Integer, power of 2) 


Parameter: 

Allowable 

Responses: 


window-choice Choice for Axial Data Window type 
0 No window (Rectangular) 

1 Tukey with 8-point taper 

2 Hanning (cosine squared) 

3 Hamming (cosine) 


Parameter: izadd 
Allowable 1 
Responses: 0 


1st coordinate (axial) Zero Padding 

Yes, zero-pad input data 

No, do not zero-pad input data 


Parameter: ihstaxt 
Allowable “f’ 
Responses: 


Starting hologram index (1st frequency bin) 
Index of first frequency bin to be processed 
from the available set (positive integer) 


Parameter: ihstep 
Allowable “m” 
Responses: 


Hologram index step interval 
Increment to the next successive frequency 
bin to be processed from the available set 
(positive integer) 
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Parameter: ihno 
Allowable “n” 
Responses: 

Parameters: rkc, alpha 
Allowable “0, 0” 
Responses: “-1, 0” 

“a, 6” 

Parameter: wname 
Allowable “name. ext” 
Responses: 


Parameter: nfilt 
Allowable 0 
Responses: 

1 

Parameter: nf jaax 
Allowable “nmai” 
Responses: 

Parameter: ifqout 
Allowable 0 
Responses: 1 
2 

3 

4 

5 


Number of holograms to process 

Total number of frequency bins to be processed 

from the available set (positive integer) 


Wavenumber domain filter parameters 
Do not apply any wavenumber filter 
Read in a frequency-dependent filter from an 
external file 

K(cutoff) radians per length, Slope 


Frequency-dependent wavenumber filter File Name 

Full path and file name of file 

which contains parameters for filter set. 

(String, < 80 characters) 


Circumferential Order Filtering flag 
Include all available circumferential orders 
(No filtering) 

Remove circumferential orders above a set maximum 


Maximum circumferential order cutoff (nfilt = 0 only) 
The maximum circumferential order included 
in reconstruction (positive integer, < ncirpts/2 


Genah Output Processing Type 
Normal Reconstruction (P, V, I) 

K-space velocity without propagation (K) 

K-space velocity with propagation (G) 

V{z,n) (NO inverse transform on 0) (N) 

P{k^,<f)) {no propagation, no transform on 0) (Q) 
V{z, 0') for orders n = n^i„,..., nmax- 
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Parameters: nmin, nmax 

Minimum and Maximum Circumferential 


orders to keep in processing 

j4iiowabie 

Lowest and highest orders to retain 

Responses: 

(positive integers; 0 < rimin < Umax < ncirpts/2 

Parameter: ichng 

Processing parameters change flag 

Allowable 0 

The reported parameters are correct: 

Responses: 

Proceed with reconstruction. 

1 

The reported parameters need to be altered. 

Parameter: alat 

Axial Lattice Spacing (ichng = 0 only) 

Allowable “a” 

Axial lattice spacing in units “units” 

Responses: 

(Positive real number) 

Parameter: drhx 

Standoff Distance (ichng = 0 only) 

Allowable “6” 

Distance, in “units”, from hologram 

Responses: 

surface to the desired reconstruction surface 


(Positive real number) 

Parameter: radius 

Reconstruction Radius (ichng = 0 only) 

Allowable “c” 

The radius of the desired reconstruction surface 

Responses: 

in units “units” (Positive real number) 

Parameter: fnamel 

Output File Name 

Allowable name. *mn 

mn = 10 X m.n = 10 x kr[cutoff) 

Responses: name.*vw 

“vw” for variable window. 


frequency-dependent wavenumber filter) 
(Note that is literal here) 

(String, < 80 characters) 
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Appendix D 

List of Routines in GEN AH 


This section contains a list of all the routines needed to run Genah . Except 
for genah itself, which appears first, they are listed in alphabetical order. 

Each description contains a short definition, any important points to be 
noted, lists of arguments, input and output parameters and routines called. 

The latter gives the file dependencies and are listed in the order in which 
they appear in the calling routine. All routines are contained in separate 
files. 

Note that only the master Bessel function routines DBESJ, DBESY and 
DBESK are listed; their dependent subroutines are not. They are all members 
of the VAXMATH package, where they are fully annotated. To repeat them here 
would waste space and detract from the intent of listing Genah routines. 
Likewise, the many implicit FORTRAN-77 functions called by Genah are not 
listed. 

program genah Routine for processing Cylindrical Nearfield Acoustical Holo¬ 
grams. 

Arguments: None 

Input Parameters; As described in Appendix C 
Output Parameters: As described in Chapter 4 

Routines Called: filread.ieee, tukey, banning, write_to_ihead, 
filwrite.ieee, wtrev, fft, windo, cang, DBESJ, DBESY, DBESK 

real function cang Function to compute the phase angle (in degrees) of a 
complex number 
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Note: This function is the analog of the FORTRAN built-in function 
cabs. An offset of 10"^^ is added to each part to obviate the 
indeterminate case tan" ^(0/0). Note that some systems do define 
this case for their C language implementations of atan2, but the 
FORTRAN-77 standard does not. 

Arguments: z 
Input Parameters: 

z (Complexes) Complex number whose phase is to be computed 
Output Parameters: 

cang (Real*4) Phase angle (in degrees) of complex number z 
Routines Called: None 

subroutine cvt_2byte_toJnt4 Routine to convert two consecutive bytes 
into an integer*4 value 

Arguments: inti, bytel, ieee 
Input Parameters: 

bytel (Byte array, 2) Consecutive bytes to be converted 
ieee (Integer*4) Number format conversion flag 

1 Convert data from DEC to IEEE number format 
0 Data is not converted (can be DEC or IEEE) 

Output Parameters: 

inti (Integer*4) Resulting integer value 
Routines Called: None 

subroutine cvt_4byte_to_real4 Routine to convert four consecutive bytes 
into a real*4 value 

Arguments: rl, bytel, ieee 
Input Parameters: 

bytel (Byte array, 4) Consecutive bytes to be converted 
ieee (Integer*4) Number format conversion flag 

1 Convert data from DEC to IEEE number format 
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0 Data is not converted (can be DEC or IEEE) 

Output Parameters: 

rl (Real*4) Resulting single precision floating point value 
Routines Called: None 

subroutine cvt_decf Jeeef Routine to convert an array of DEC single pre¬ 
cision floating point numbers to IEEE floating point. 

Arguments: dec_fp, ieee_fp, num 
Input Parameters: 

dec_fp (Real*4, num) Array of DEC floating-point numbers 
num (Integer*4) Number of values to be converted 
Output Parameters: 

ieee_f p (Real*4, num) Array of converted, word-swapped IEEE 
formatted numbers 

Routines Called: None 

subroutine cvt Jnt4_to_2byte Routine to convert an integer*4 value into 
two consecutive bytes 

Arguments: inti, bytel, ieee 
Input Parameters: 

inti (Integer*4) Integer value to be converted 
ieee (Integer*4) Number format conversion flag 

1 Convert data from DEC to IEEE number format 
0 Data is not converted (can be DEC or IEEE) 

Output Parameters: 

bytel (Byte array, 2) Consecutive bytes converted from integer 
Routines Called: None 

subroutine cvt_real4_to_4byte Routine to convert a real*4 value into four 
consecutive bytes 

Arguments: rl, bytel, ieee 
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Input Parameters: 

rl (Real*4) Single precision floating point value to be converted 
ieee (Integer*4) Number format conversion flag 

1 Convert data from DEC to IEEE number format 
0 Data is not converted (can be DEC or IEEE) 

Output Parameters: 

bytel (Byte array, 4) Consecutive bytes converted from floating 
point value 

Routines Called: None 

subroutine DBESJ Routine to compute an n member sequence of the dou¬ 
ble precision Bessel Function of the First Kind, Ja{x) for non-negative 
a and double precision argument. 

Note: VAXMATH routine 

Arguments: darg, alpha, maxn, djarray, nzero 
Input Parameters: 

darg (Real*8) Double precision argument, x > 0 
alpha (Real*8) Double precision order of first member in sequence, 
a > 0 

majcn (Integer*4) Number of members in sequence, n > 1 
Output Parameters: 

dj array (Real*8) Double precision vector result; the first meixn 
components contain values for Ja+k-i(^), A; = 1,... ,maxn 
nzero (Integer*4) Number of components of dj array set to zero 
due to underflow. 

Routines Called: DASYJY, DJAIRY, DLNGAM, DIMACH, IlMACH, XERROR, 
all from VAXMATH 

subroutine DBESY Routine to compute an n member sequence of the 
double precision Bessel Function of the Second Kind, Ya{x) for non¬ 
negative a and double precision argument x. 

Note: VAXMATH routine 
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Arguments: daxg, alpha, maxn, dyarray 
Input Parameters: 

daxg (Real*8) Double precision argument, a: > 0 
alpha (Real*8) Double precision order of first member in sequence, 
a > 0 

maxn (Integer*4) Number of members in sequence, n > 1 
Output Parameters: 

dyarray (Real*8) Double precision vector result; the first maxn 
components contain values for yQ+jt-i(x), A: = 1,... ,maxn 

Routines Called: DASYJY, DBESYO, DBESYl, DBSYNU, DYAIRY, DIMACH, 
IlMACH, XERROR, all from VAXMATH 

subroutine DBESK Routine to compute an n member sequence of the 
double precision Modified Bessel Function of the Second kind, Ka{x) 
for non-negative a and double precision argument x. 

Note: VAXMATH routine 

Arguments: darg, alpha, if lag, maxn, dkarray, nzero 
Input Parameters: 

daxg (Real*8) Double precision argument, x > 0 

alpha (Real*8) Double precision order of first member in sequence, 

Q > 0 

iflag (Integer*4) Scaling option flag; iflag = 1 returns Ka{x), 
iflag = 2 returns e^Ka{x) 

maxn (Integer*4) Number of members in sequence, n > 1 
Output Parameters: 

djarray (Real*8) Double precision vector result; the first maxn 
components contain values for Ka+k-iix), k = 1,... ,maxn 
nzero (Integer*4) Number of components of djaxray set to zero 
due to underflow when iflag = 1 

Routines Called: DASYIK, DBESKO, DBESKl, DBSKNU, DBSKOE, DBSKIE, 
DIMACH, IlMACH, XERROR, all from VAXMATH 
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subroutine fFt Routine to compute the Fast Fourier Transform of a complex 
vector 

Note: Transform is computed in place 
Arguments: isw, n, 12n, ibr, i2, w, z 
Input Parameters: 

isw (Integer*4) Forward/Inverse transform switch 

0 

1 

n (Integer*4) Length of transform vector 
12n (Integer*4) /op 2 (n) (number of powers of 2) 
ibr (Integer*4, n) Bit reversal index vector 
i2 (Integer*4, 12n) Power of 2 index vector 
w (Complex*8, n/2) Exponential weighting factors 
z (Complexes, n) Vector to be transformed 
Output Parameters: 

z (Complex*8, n) Coefficients of transformed vector 
Routines Called: None 

subroutine filreadJeee Routine to open, read and close Genah format 
data files. 

Arguments: iopen, lu, ieee, freq, filename, ihol, Idz, zc 
Input Parameters: 

iopen (Integer*4) Routine control fiag 

1 Open the file and read header, put header information into 
the common block. 

0 Read hologram number ihol and load data into zc 
-1 Close the input file 

lu (Integer*4) Logical unit for file to be opened/read 
ieee (Integer*4) Number format conversion flag 

1 Convert data from DEC to IEEE number format 
0 Data is not converted (can be DEC or IEEE) 
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Idz (Integer*4) Leading dimension of data array zc 
ihol (Integer*4) Hologram number corresponding to current record(s) 
Output parameters: 

freq (Real*4) Frequency corresponding to ihol 
filename (Character*80) Character string containing file name, 
passed back after OPEN 

zc (Complex*8) Array of data values for bin ihol of file filename 
Routines Called: read_from_ihead, cvt_decf_ieeef 

subroutine filwriteJeee Routine to open, write and close Genah format 
data files. 

Arguments: iopen, lu, ihead, num_recs, filename, ihol, Idz, 
zc 

Input Parameters: 

iopen (Integer*4) Routine control flag 

1 Open the output file and write the header block. 

0 Write hologram number ihol from array zc 
-1 Close the output file 

lu (Integer*4) Logical unit for file to be opened/written 
ihead (Byte) 512-byte array with encoded header information 
nimi_recs (Integer*4) Number of records to write, equivalent to 
number of scan lines (ncirpts) 

filename (Character*80) Character string containing file name, 
passed back after OPEN 

ihol (Integer*4) Hologram number corresponding to current record(s) 
Idz (Integer*4) Leading dimension of data array zc 
zc (Complex*8) Array of data values for bin ihol of file filename 
Routines Called: None 

subroutine banning Routine to compute the coefficients for a Hanning 
window for an even sequence 

Note: The Hanning window is one half-cycle of 
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Arguments: window, npts, order 

Input Parameters: 

npts (Integer*4) Length of window vector 
order (Integer*4) Power to which Cosine is raised 

Output Parameters: 

window (Real*4, npts) Coefficients of Hanning window of length 
npts 

Routines Called: None 

subroutine readJFromJhead Routine to read the currently open file’s header 

parameters from the byte array ihead into the header variables 

Note: Call this routine after ihead is loaded from record 1 of a broad¬ 
band holography data file. 

Arguments: ihead, ieee 

Input Parameters: 

ieee (Integer*4) Number format conversion flag 

1 Convert data from DEC to IEEE number format 
0 Data is not converted (can be DEC or IEEE) 
ihead (Byte array, 512) Raw bytes from record 1 of data file, 
declared in include file header_ieee.f 

Output P^lrameters: 

header_ieee_chax Common block containing all header param¬ 
eters which are character strings, declared in include file 
header_ieee.f 

header.ieee Common block containing all header parameters which 
are not character strings, declared in include file header_ieee. f 

Routines Called: cvt_2byte_to_int4, cvt_4byte_to_real4 

subroutine tukey Routine to compute the coefficients for a flat-top Tukey 

window 

Note: This version is hardwired to produce an 8 point taper. 
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Arguments: window, npts 
Input Parameters: 

npts (Integer*4) Length of window vector 
Output Parameters: 

window (Real*4, npts) Coefficients of 8 point taper Tukey window 
of length npts 

Routines Called: None 

real function windo Routine to compute the modulus of the two dimen¬ 
sional wavenumber domain filter 

Arguments: rkc, efold, xk, yk 
Input Parameters: 

rkc (Real*4) Cutoff wavenumber (radial) 

efold (Real*4) Exponential folding rate (steepness of the filter) 

xk (Real*4) First coordinate wavenumber {kx or k^) 

yk (Real*4) Second coordinate wavenumber {ky or n/a) 

Output Parameters: 

windo (Real*4) Modulus of filter at kj. = yjk^ + k^ 

subroutine write.to Jhead Routine to write the current header parame¬ 
ters from the header variables into the byte array ihead 

Note: Call this routine before ihead is written to record 1 of a broad¬ 
band holography data file. 

Arguments: ihead, ieee 

Input Parameters: 

ieee (Integer*4) Number format conversion flag 

1 Convert data from DEC to IEEE number format 
0 Data is not converted (can be DEC or IEEE) 
header_ieee_char Common block containing all header param¬ 
eters which are character strings, declared in include file 
header.ieee.f 
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header.ieee Common block containing all header parameters which 
are not character strings, declared in include file header.ieee. f 

Output Parameters: 

ihead (Byte array, 512) Raw bytes for record 1 of data file, de¬ 
clared in include file header.ieee.f 

Routines Called: cvt_int4_to_2byte, cvt_real4_to_4byte 

subroutine wtrev Routine to compute Fast Fourier Transform weighting 
factors and bit reversal index 

Arguments: isw, n, 12n, ibr, i2, w 
Input Parameters: 

isw (Integer*4) Forward/Inverse transform switch 

0 

1 

n (Integer*4) Length of transform vector 
12n (Integer*4) log 2 {n) (number of powers of 2) 

Output Parameters: 

ibr (Integer*4, n) Bit reversal index vector 
i2 (Integer*4, 12n) Power of 2 index vector 
w (Complexes, n/2) Exponential weighting factors 
Routines Called: None 
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Appendix E 

Wavenumber Filters in 
GENAH 


This section is intended to provide a concise description of the matched 
exponential filter used by Genah . 

Two parameters control the filter. The first, rkc, is the cutoff wave- 
number, kr(cutoff)- This is the radius of the circle on the two dimensional 
wavenumber domain where the filter reaches a value of 1/2 (or the “3 dB 
down point”, as filter designers would call it). The “radial” (actually a 
trace) wavenumber for any frequency is defined as follows: 

27r 
At * 

i,Ak^ ; -N,/2 < i, < Nj2 

n 

, ^ n <C. Tlmax/'^ 

ro _ 

where is the number of points in the axial direction, is the lattice 
spacing in the axial direction, n is the circumferential order, Umax is the 
maximum circumferential order and ro is the radius of the reconstruction 
surface. The filter magnitude for a given radial wavenumber is then defined 
as 


Ak, = 
= 

k^ = 

kj. = 


/(fcr) — 1 2 ' k,. K. kj.(^^toff) 
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kf. ^ fcr(cuto//) 


f! 

2 ’ 

where 

^ = - (1 - 7-) - 

^ ^r(cutof f) 

The second parameter controls the rate at which the filter transitions 
from 1 to 0 as a function of radial wavenumber. It is known as the steepness 
parameter or the exponential folding rate. As with any filter response, the 
steeper the roll-off, the higher the sidelobes or ringing in the inverse transform 
domain. Figure E.l shows filters for a range of steepness parameters, a. 
Typically, a = 0.05 is a good choice for most reconstructions. 

When trying to determine values for the cutoff wavenumber, use of Pro¬ 
cessing Options 1 or 2 can be valuable (see Chapter 4). Plots of the wavenum¬ 
ber domain response will reveal a relatively circular region of low response. 
Inside this region lies the data from the shell itself. Outside this region 
lies high wavenumber noise, amplified by the propagator. Placing the cutoff 
wavenumber at the radius of this circular region is most effective. Since the 
valid data region often grows radially with increasing frequency, it is nec¬ 
essary to use the “Variable Window” option, specifying new filters as the 
frequency increases. 
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Wavenumber Filter Amplitude 



Figure E.l; Plot of a set of wavenumber domain filter responses for 
values of the steepness parameter, a. 
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